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ABSTRACT 

We present a series of models for the plasma properties along open magnetic flux tubes rooted in solar coronal 
holes, streamers, and active regions. These models represent the first self-consistent solutions that combine: 
(1) chromospheric heating driven by an empirically guided acoustic wave spectrum, (2) coronal heating from 
Alfven waves that have been partially reflected, then damped by anisotropic turbulent cascade, and (3) solar 
wind acceleration from gradients of gas pressure, acoustic wave pressure, and Alfven wave pressure. The only 
input parameters are the photospheric lower boundary conditions for the waves and the radial dependence of 
the background magnetic field along the flux tube. We have not included multifluid or collisionless effects 
(e.g., preferential ion heating) which are not yet fully understood. For a single choice for the photospheric 
wave properties, our models produce a realistic range of slow and fast solar wind conditions by varying only 
the coronal magnetic field. Specifically, a two-dimensional model of coronal holes and streamers at solar 
minimum reproduces the latitudinal bifurcation of slow and fast streams seen by Ulysses. The radial gradient 
of the Alfven speed affects where the waves are reflected and damped, and thus whether energy is deposited 
below or above the Parker critical point. As predicted by earlier studies, a larger coronal "expansion factor" 
gives rise to a slower and denser wind, higher temperature at the coronal base, less intense Alfven waves at 1 
AU, and correlative trends for commonly measured ratios of ion charge states and FlP-sensitive abundances that 
are in general agreement with observations. These models offer supporting evidence for the idea that coronal 
heating and solar wind acceleration (in open magnetic flux tubes) can occur as a result of wave dissipation and 
turbulent cascade. 

Subject headings: MHD — solar wind — Sun: atmospheric motions — Sun: corona — turbulence — waves 



1. INTRODUCTION 

One of the most persistent problems in solar physics has 
been the unambiguous identification of the mechanisms that 
heat the Sun's corona and accelerate the solar wind. Many 
processes have been proposed for converting some fraction 
of the mechanical energy in subphotospheric convective mo- 
tions to heat, but it has proved very difficult to make distin- 
guishing comparisons between the predictions of these com- 
peting ideas and specific observations. We are entering an 
era, though, where both the models and the measurements are 
improving to the point of soon being able to eliminate many 
of the candidate theories. For example, it seems increasingly 
clear that closed loops in the low corona are heated by small- 
scale, intermittent magnetic reconnection that is driven by the 
continual stressing of their magnetic footpoints (see recent re- 
views by Longcope 2004; Gudiksen 2005; Aschwanden 2006; 
Klimchuk 2006). 

In this paper we model the coronal heating along open field 
lines that reach into interplanetary space. We construct a self- 
consistent model of the photosphere, chromosphere, corona, 
and solar wind that is driven mainly by magnetohydrody- 
namic (MHD) turbulence and is free of arbitrary "heating 
functions" that have been used in the past. This work contin- 
ues earlier studies (Cranmer & van Ballegooijen 2005; Cran- 
mer 2005a) which used prescribed empirical descriptions of 
the density and flow velocity along an open flux tube in or- 
der to compute the rates of turbulent heating and acceleration. 
The main goal of this paper is to provide a possible explana- 
tion for the origin and properties of fast and slow solar wind 
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streams. Another goal, though, is to illustrate how the partic- 
ular description of MHD turbulence can be applied to more 
advanced modeling efforts — specifically those that attempt to 
reproduce the full three-dimensional and time-dependent na- 
ture of the corona and heliosphere (e.g., Riley et al. 2001; 
Roussev et al. 2003; Toth et al. 2005; Usmanov & Goldstein 
2006). 

Coronal heating and solar wind acceleration have been 
known to be closely linked since the initial contributions of 
Parker (1958). However, nearly all subsequent theoretical at- 
tempts to model both processes together have made limiting 
assumptions about either the heating (i.e., ad hoc energy in- 
put rates) or the acceleration (i.e., a prescribed mass flux). It 
has been realized over the past few decades that it is key to 
resolve the full chromosphere-corona transition region in or- 
der to produce a model with an internally consistent radial 
dependence of pressure, density, and flow speed (see Ham- 
mer 1982; Hansteen & Leer 1995; Hansteen et al. 1997; Lie- 
Svendsen & Esser 2005). To our knowledge, the only solar 
wind models — other than the ones presented in this paper — 
that contain both a first-principles approach to coronal heat- 
ing and a self-consistent chromosphere, transition region, and 
mass flux are those of Suzuki & Inutsuka (2005, 2006). 1 

The study of MHD turbulence as a potential source of heat- 
ing for the solar wind goes back to Coleman (1968) and 
Jokipii & Davis (1969). Hollweg (1986) extended these ideas 

1 Other models can be included if some of the above conditions are relaxed. 
For example, Hu et al. (2000) and Li (2002, 2003) considered wave-driven 
coronal heating with a lower boundary within the transition region (i.e., tem- 
peratures ranging between 6 X 10 4 and 8 X 10 5 K). See § 2 below for a com- 
paidgaiSiffliiFiaialaKiBtipliysical assumptions and numerical approaches. 
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down into the corona and laid the foundations for the cascade- 
driven heating rates used by Isenberg (1990), Li et al. (1999), 
Matthaeus et al. (1999), Dmitruk et al. (2001, 2002), and this 
paper. These ideas are highly complementary to theoretical 
models of wave dissipation via ion cyclotron resonance (e.g., 
Hollweg & Isenberg 2002; Cranmer 2000, 2001, 2002) that 
have been invoked to explain the observed preferential heat- 
ing of heavy ions in the extended corona (Kohl et al. 1997, 
1998, 2006). In the cascade paradigm, the energy in large- 
scale Alfvenic fluctuations must eventually be dissipated in 
small-scale (short wavelength or high frequency) collisionless 
kinetic modes. The results from this paper can thus be used 
as initial conditions for detailed models of the anisotropic tur- 
bulent cascade and minor ion heating. 

The remainder of this paper is organized as follows. In § 2 
we present an overview of the general principles involved in 
our modeling effort (i.e., why certain physical processes are 
included or excluded). § 3 gives the conservation equations 
that are solved and the adopted prescriptions for basic ingre- 
dients such as radiative heating/cooling and conduction. The 
detailed treatment of waves is described in the following three 
sections: acoustic waves and shock heating in § 4, Alfven 
waves and turbulent cascade in § 5, and the ponderomotive 
wave pressure acceleration due to both types of fluctuations 
in § 6. The numerical methods used to solve the equations 
and the newly developed computer code, called ZEPHYR, are 
described in § 7. We then present a comprehensive set of solu- 
tions for the solar wind emerging from coronal holes, helmet 
streamers, and active regions (§ 8). A summary of the major 
results of this paper, together with a discussion of the impli- 
cations for understanding the winds of other stars, is given in 
§ 9. 

2. MODEL PHILOSOPHY 

We consider the one-dimensional variation of plasma pa- 
rameters along a radially pointed magnetic flux tube rooted in 
the solar photosphere and open to interplanetary space. Some 
effects of the multi-dimensional magnetic field in the photo- 
sphere and chromosphere are included as explicit superradial 
expansion of the flux tube (see also Cranmer & van Balle- 
gooijen 2005), but otherwise we ignore the effects of neigh- 
boring closed flux tubes. We assume that most of the plasma 
that eventually becomes the time-steady solar wind originates 
in small (100 km sized) intergranular flux concentrations that 
are concentrated most densely in the supergranular network. 

The models we construct are time-independent solutions to 
the hydrodynamic conservation equations. This may be a re- 
strictive simplification since all observed layers of the solar at- 
mosphere are intrinsically dynamic over a wide range of time 
scales. It is evident, for example, that the transition region is 
"strongly nonuniform and magnetically structured" (Marsch 
et al. 2006) to the extent that one-dimensional layered mod- 
els may miss some important physical effects at those heights. 
However, the corona above the transition region seems to ex- 
ist in its hot (r ~ 10 6 K) state all the time. Also, in situ mea- 
surements of the solar wind plasma often can be interpreted 
clearly as a superposition of a quasi-steady supersonic flow 
(which depends mainly on the coronal source region of the 
streams being measured or on stream-stream interactions) and 
a rapidly varying turbulent or wavelike component. We thus 
separate these two scales and model the fluctuations using en- 
ergy conservation equations derived from linear wave theory. 

It is especially beneficial to treat waves and turbulent mo- 
tions statistically — rather than follow the oscillations explic- 



itly in time — when the dynamically important ranges of wave- 
length and period span many orders of magnitude. Alfven 
waves in the solar wind are measured to have periods from 
seconds to days (e.g., Goldstein et al. 1995; Tu & Marsch 
1995) and there is evidence that the waves that dominate the 
dissipation in the collisionless extended corona have periods 
even below 10~ 3 s (McKenzie et al. 1995; Cranmer et al. 
1999). If this full range of scales had to be resolved in order 
to track the motions in a time-dependent numerical model, 
it would be extremely difficult to simulate the macroscopic 
plasma properties simultaneously from the photosphere to the 
heliosphere. 

There are, however, some clear disadvantages in mod- 
eling fluctuations as statistical "wave energy fluxes" rather 
than as explicit variations. Any nonlinear processes, such 
as shock steepening, turbulent cascade, mode conversion, 
or ponderomotive forces, must be inserted beforehand as 
source or sink terms in the time-steady conservation equa- 
tions. Time-dependent MHD models have the benefit of pro- 
ducing such effects naturally (e.g., Ofman & Davila 1998; 
Ofman 2005; Bogdan et al. 2002, 2003; Suzuki & Inutsuka 
2005, 2006). However, by resolving only certain temporal 
and spatial scales, these models are limited in ways that a 
statistical treatment is not. For example, the relative impor- 
tance of shock dissipation in the time-dependent models of 
Suzuki & Inutsuka (2005, 2006) may be an artifact of ei- 
ther the inability to resolve small enough scales (which could 
drive processes such as perpendicular cascade, Landau damp- 
ing, and collisionless particle energization at shocks) or the 
waveguide-like trapping of fluctuations along the model flux 
tube. It is important to note, though, that neither extreme — 
i.e., neither coarse time-dependent simulations nor statistical 
(and more approximate) time-independent models — is ideal. 
There may be value in seeking some kind of hybrid method- 
ology between these two approaches. 2 

As indicated above, the only source of heating for the chro- 
mosphere and corona that we consider is the dissipation of 
waves and turbulent motions. For open flux tubes, a sub- 
stantial fraction of the energy appears to be deposited at 
large heights in the wind's acceleration region — i.e., at spatial 
scales much larger than the sizes of low-lying closed loops 
in the quiet Sun and active regions. This demands the en- 
ergy be propagated for some distance, presumably by waves 
or turbulent eddies, before it is dissipated. However, the gen- 
eral phenomenology of turbulence is probably not limited to 
the open-field regions. Concepts from turbulence theory have 
been applied to the full range of time scales for closed-field 
coronal energy input as well, from the most rapid (AC) wave- 
like oscillations to the slowest (DC) quasistatic stresses on 
magnetic footpoints (see van Ballegooijen 1986; Hendrix & 
van Hoven 1996; Milano et al. 1997; Gomez et al. 2000; Chae 
et al. 2002). Indeed, some recent simulations of intermittent 
turbulent heating in closed loops have been interpreted using 
similar cascade rate expressions as we use in this paper (Rap- 
pazzo et al. 2007). 

The models presented below include the effects of both 
Alfven waves and acoustic waves on the mean flow. We as- 
sume implicitly that all waves propagate parallel to the ra- 
dially oriented flux tube, but this is not an essential feature. 
In the magnetically dominated corona (where /?< 1; /3 be- 

2 Recent advances in modeling MHD turbulence in coronal loops with a 
so-called "shell model" in wavenumber space may be pointing the way (e.g., 
Giuliani & Carbone 1998; Nigro et al. 2004). 
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ing the ratio of gas pressure to magnetic pressure) the mod- 
eled Alfven waves can be considered essentially to be the 
sum of true Alfvenic wave power and fast-mode MHD wave 
power. The acoustic waves can be considered equivalent to 
slow-mode MHD waves. As a starting point, we neglect all 
nonlinear couplings between the Alfven and acoustic wave 
modes. By ignoring the enhanced reflection and dissipation 
that may arise because of these couplings (e.g., Suzuki & In- 
utsuka 2005, 2006) we essentially underestimate the heat de- 
posited in the extended corona and heliosphere. The fact that 
sufficient energy nonetheless exists to heat the corona and ac- 
celerate the wind (mainly from the Alfven waves) appears to 
deemphasize the need for such mode couplings. 

It should be noted that the attention given to modeling the 
steepening and shock dissipation of the acoustic waves (§ 4) 
does not seem to have much of a "payoff" in the resulting 
properties of the corona and solar wind. The fully ionized 
outer atmosphere is extremely insensitive to the magnitude 
(or sometimes even the presence) of acoustic wave power 
that comes up from below the photosphere (§ 8.2). Accord- 
ingly, many prior studies of the coupled chromosphere and 
corona used a much simpler prescription for maintaining the 
chromosphere. We believe it is important, though, to treat 
both the chromospheric and coronal heating at a comparable 
level of first-principles modeling. We also anticipate that the 
ZEPHYR code developed here will be applied to the simula- 
tions of winds of other late-type stars that may have acousti- 
cally heated coronae (e.g., Mullan & Cheng 1993; Schrijver 
1995), and thus the acoustic waves should be treated as real- 
istically as possible. 

One final limitation of the models presented below is that 
the energy conservation is treated in a one-fluid manner. The 
protons, electrons, and heavy ions are modeled as having a 
common flow speed u and temperature T, with microscopic 
velocity distributions that are simple isotropic Maxwellians. 
This is an extreme simplification of reality, since it has been 
known since the 1960s that the in situ solar wind exhibits 
significant departures from a thermalized equilibrium (see re- 
views by Hundhausen 1972; Feldman & Marsch 1997). Many 
of these effects persist down into the extended corona as well 
(e.g., Kohl et al. 2006). Therefore, theoretical models have 
long included a range of attempts to deal with these features. 
Fluid-based models have been extended to solve separate con- 
servation equations for each particle species, and they have 
been reconstructed in various ways based on different param- 
eterizations for the anisotropic velocity distributions. Several 
purely kinetic models have also been attempted (see above- 
cited reviews, also Cranmer 2002; Hollweg & Isenberg 2002; 
Marsch 2005). 

Despite the potentially clumsy "averaging" over real kinetic 
effects, we believe the one-fluid approach is the most consis- 
tent with our present state of knowledge about the primary 
source of energy deposition: MHD turbulent cascade. There 
are a number of competing ideas in the literature regarding 
how the cascade proceeds to its smallest kinetic scales and 
how either linear or nonlinear damping transfers wave energy 
to the particles (for a recent summary see § 5.2.4 of Kohl et al. 
2006). We thus do not yet know, from first principles, how to 
partition the cascaded wave energy between particle species 
(T e y^T p ^ 7i on ) and between various directions in microscopic 
velocity space (7]| =/ 7j_). Performing such partitioning at 
the present time would essentially add new free parameters 
into the model. Our adopted rate of Alfvenic coronal heating 
(eq. l47l ) thus deals only with the total energy flux that cas- 



cades from large to small scales and not the specific means of 
dissipation once the energy reaches the small scales. 

3. BASIC PHYSICS 

The equations governing the expansion of a time-steady 
stellar wind are derivable by taking successive velocity mo- 
ments of the Boltzmann equation, together with some as- 
sumption about the shape of the velocity distribution function 
in order to close the otherwise infinite chain of moment equa- 
tions (see, e.g., Braginskii 1965; Collins 1989; Marsch 2005). 
In this section we describe the conservation equations used 
in the models (§ 3.1) and the adopted prescriptions for heat 
transport due to radiation (§ 3.2) and conduction (§ 3.3). 

3.1. Conservation Equations 

We consider the flow of a pure hydrogen plasma along a 
radially oriented magnetic field. The goal is to solve for the 
time-steady radial dependence of the mass density p, the bulk 
flow speed u, and the Maxwellian temperature T. The dis- 
tance along the magnetic flux tube is denoted either as r, mea- 
sured from Sun-center, or z, measured from the lower bound- 
ary of the model in the solar atmosphere (essentially the pho- 
tosphere). For completeness, the equations below contain the 
dependence on time t , though these terms are set to zero in the 
time-independent solutions that we describe below. 

The equation of mass conservation is 



(i) 



where A is the cross-sectional area of the one-dimensional 
flux tube along which the wind flows. Magnetic flux conser- 
vation demands that the product BqA is constant along the flux 
tube, where Bo is the field strength that we specify explicitly 
(see § 8). A time-steady one-dimensional flow thus constrains 
the product pu to be proportional to Bo- 
The equation of momentum conservation is 



du du 1 dP _ GM* 
dt dr p dr r 1 



+ D 



(2) 



where P is the gas pressure, G is the Newtonian gravitation 
constant, and is the mass of the star. The mass of the 
plasma in the modeled stellar wind region is assumed to be 
negligible from a gravitational standpoint. Also, D is the bulk 
acceleration on the plasma due to wave pressure; i.e., the 
nondissipative net ponderomotive force due to the propaga- 
tion of waves through an inhomogeneous medium (§ 6). 
The equation of internal energy conservation is 

dE BE [E + P\ d , „ „ „ „ _ 

— + U — + I -j- J — (l*A) = Qml + Qeond + QA + Qs (3) 

where E is the internal energy density and the terms on the 
right-hand side are volumetric heating/cooling rates due to 
radiation, conduction, Alfven wave damping, and acoustic 
(sound) wave damping. The terms on the left-hand side that 
depend on u are responsible for enthalpy transport and adia- 
batic cooling in the accelerating wind. The coupling of the 
above equations requires additional constitutive relations to 
be specified: 

P = n tol k B T , (4) 



■+n p In , 



7-1 

«tot = n H + n e = (n Q + n p ) + n e , 



(5) 
(6) 
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FIG. 1. — Temperature dependence of the hydrogen ionization state x 
(dashed line), the corresponding neutral hydrogen fraction 1 — x (dotted line), 
and the optically thin radiative loss function A(T) in units of 10~ 22 erg cm 3 
s _1 (solid line). 

where k% is Boltzmann's constant and a monatomic ratio of 
specific heats 7 = 5/3 is used. The total particle number den- 
sity n tot is given by the sum of the hydrogen number density 
«h and the electron density n e , and the number densities of 
neutral hydrogen and protons are denoted no and n p , respec- 
tively. For the pure hydrogen plasma assumed in this paper, 
n p = n e . The mass density is given by p = mftnn + m e n e , al- 
though the second term is safely considered to be negligible. 

Note that equation <(5j becomes the ideal gas relation E = 
3P/2 = pc v T in the purely neutral limit, where c v is the spe- 
cific heat at constant volume. Ionization is taken into account 
by the internal energy's dependence on In, the ionization po- 
tential of hydrogen from the ground level (13.6 eV). We use 
the same convention in the definition of E as Ulmschneider & 
Muchmore (1986) and Mullan & Cheng (1993). An alternate 
definition of the internal energy is possible, though, which 
is essentially given by E-nnln- This version reduces to the 
above ideal-gas relation in the limit of a fully ionized plasma 
(e.g., McClymont & Canfield 1983; Fontenla et al. 1990). Ei- 
ther definition is consistent with the combined equations of 
mass and energy conservation given above. 

We adopted a relatively simple prescription to compute the 
ionization fraction x = n p /nn- This quantity is parameterized 
as a tabulated function of T only, where we have used the 
ionization balance from a recent semi-empirical model of the 
solar photosphere, chromosphere, and transition region (E. 
Avrett 2005, private communication; see also Fontenla et al. 
1993, 2002, 2006). The tabulated model we use is a modi- 
fied version of the quiet-Sun model C from the above papers, 
and it is the same one used by Cranmer & van Ballegooijen 
(2005) to model the magnetic structure of the chromospheric 
network. For convenience we call this model FAL-C'. 

Figure 1 shows the adopted ionization fraction x, as well as 
the corresponding neutral fraction 1 —x (in order to see small 
departures from total ionization at high temperatures) and the 
radiative cooling function that was derived in part from the 
FAL-C' model (see § 3.2). For temperatures below the min- 
imum tabulated value of about 4500 K, we prevent x from 
decreasing below a minimum value of 3 x 10 -4 , because pho- 



toionization is expected to always keep some small fraction 
of metals ionized at and above the photosphere. 

We tested the ZEPHYR code with a more self-consistent 
model of the hydrogen ionization balance. The results were 
extremely similar, though, to those using the tabulated frac- 
tion shown in Figure 1. This model consisted of a three- 
level hydrogen atom, where the n = 1 and n = 2 levels were 
assumed to remain in relative local thermodynamic equilib- 
rium (LTE) and the full rate equation between n = 2 and the 
continuum was solved iteratively with collisional and radia- 
tive terms from Hartmann & MacGregor (1980), Vriens & 
Smeets (1980), Ferland et al. (1992), and Ferguson & Fer- 
land (1997). Because the tabulated version was much faster 
in terms of computation time, though, we decided to use it in 
the models shown below. It will be important to include this 
kind of self-consistent ionization balance when adapting this 
method to the winds of other stars (e.g., Natta et al. 1988). 

3.2. Radiative Heating and Cooling 

A complete treatment of the non-LTE transfer of energy be- 
tween radiation and matter in a partially ionized plasma is 
beyond the scope of this paper. Detailed computational ef- 
forts to model chromospheric and coronal radiative transfer 
effects (e.g., Avrett & Loeser 1992; Carlsson & Stein 1997, 
2002; Rammacher et al. 2005) are important for reproduc- 
ing the spectrum or studying time-dependent or multidimen- 
sional dynamics, but for our purposes a simpler treatment is 
warranted. We use a similar general "bridging" approach as 
Mullan & Cheng (1993) to combine different limiting cases 
in the optically thick lower atmosphere and the optically thin 
upper atmosphere. 

The adopted radiative heating/cooling rate is given by 

Grad = e" TR/TO e thl n + (l -e" TR/T °)ethick (7) 

where t r is the Rosseland mean optical depth and tq = 0.1 
is a constant that defines where the rate is dominated by the 
optically thick or thin limits (see also Ludwig et al. 1994). 
For a spherical stellar atmosphere, we use a definition for the 
Rosseland optical depth, 

dr R = -K R p dr ' 1 ( 8 ) 

that contains the geometrical correction factor suggested by 
Lucy (1971, 1976); it is unimportant in the case of the Sun — 
where r w in regions where t r w 1 — but we include it for 
later use in the extended atmospheres of evolved stars. The 
Rosseland mean opacity k r (in cm 2 g" 1 ) is interpolated as 
a function of temperature and pressure from the table of Ku- 
rucz (1992). For a given one-dimensional model, we integrate 
downward from the upper boundary far out in the supersonic 
wind (which has an assumed optical depth of zero) to compute 
r R (r). 

In the optically thick photosphere and lower chromosphere 
we assume that the heating and cooling is dominated by con- 
tinuum photons emitted and absorbed in LTE, with 

Cthick = 47rp / K v (J v -S v )dv = 4npK R (J-S) (9) 

(e.g., Mihalas 1978; Vogler et al. 2004). Above, S is the 
frequency-integrated source function, assumed here in LTE 
to equal the local integrated Planck function, S = B = a R T 4 / ir, 
with or being the Stefan-Boltzmann constant. The frequency- 
integrated mean intensity J is given by the gray atmosphere 
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dependence on optical depth r R , 

/(tr) = ^^R^eff [T R + 2q(T R )W(r)] , (10) 

where T e ff = 5800 K is the solar effective temperature and q(j) 
is the Hopf (1930, 1932) function, for which we use the fol- 
lowing fit (accurate to better than 0.2%), 



q(r) = 0.710- 



0.133 



(1+0.15-r - 73 ) 



0.73U7.4 



(ID 



Also, W(r) = [l-(l-fl 2 /r 2 ) 1/2 ]/2 is the spherical dilution 
factor suggested for use in this context by Chandrasekhar 
(1934) and Lucy (1971, 1976). Qthick is positive — giving net 
heating — when T is less than the radiative equilibrium tem- 
perature (r ra d) given by J = S, and it is negative — with net 
cooling — when T > T^. 

In the optically thin upper chromosphere and corona, we 
use a modified temperature-dependent radiative cooling func- 
tion that has been computed for photon losses due to a large 
collection of spectral lines and continuum processes, 



awn = -» e » H A(r) [ l-^f 



(12) 



(see also Cox & Tucker 1969; Anderson & Athay 1989b; 
Schmutzler & Tscharnuter 1993). The radiative loss function 
A(r) is shown in Figure 1, and it has been assembled from 
three sources. (1) For hydrogen and helium above 10 4 K, we 
used the output from the PANDORA radiative transfer code 
which produced the FAL-C' semi-empirical model discussed 
above. (2) For other elements above 10 4 K, we used a tabu- 
lated radiative loss function from version 4.2 of the CHIANTI 
atomic database (Young et al. 2003) using a traditional solar 
abundance mixture (Grevesse & Sauval 1998) and collisional 
ionization balance (Mazzotta et al. 1998). (3) For temper- 
atures below 10 4 K, we used the fitting function for partially 
ionized hydrogen given by Scholz & Waters (1991). For these 
low temperatures we also added a constant lower-limit value 
of 10~ 34 erg cm 3 s" 1 to A(T) to account for photoionized met- 
als, molecules, and dust, which may be important contributors 
to radiative cooling at very low temperatures in late-type stel- 
lar atmospheres (e.g., Schirrmacher et al. 2003). 

In order to ensure that the optically thin parts of the atmo- 
sphere would smoothly approach radiative equilibrium (in the 
absence of nonradiative heating) in the same way as in the op- 
tically thick parts of the atmosphere, we multiplied the stan- 
dard cooling function by the term in parentheses in equation 
(TT2l) . For temperatures above ~ 2Z[- ac i this correction factor 
rapidly approaches unity. 

Note from Figure 1 that there is no isolated "Lyman alpha 
peak" at temperatures of 1-2 xlO 4 K. For solar atmosphere 
models, this peak seems to be a spurious feature that appears 
in optically thin radiative loss curves computed without a ra- 
diation field. The FAL-C' model used here contained a full 
non-LTE treatment of hydrogen as well as ambipolar diffusion 
that also provides additional smearing of discrete features re- 
lated to the strong H I Lya transition (see also Kuin & Poland 
1991;Fontenla et al. 2002). 

With regard to the extension to other stars, the computation 
of the radiative cooling rate is in a similar situation as the ion- 
ization fraction discussed in § 3.1. Our use of tabulated quan- 
tities from FAL-C' seems to be reasonable for modeling the 
solar atmosphere, but it will eventually need to be replaced 
with a more self-consistent procedure. Specifically, evolved 



low-gravity stars with high mass loss rates are expected to 
have more optically thick chromospheres. Hartmann & Mac- 
Gregor (1980), Canfield & Ricchiazzi (1980), and Mullan & 
Cheng (1993) computed approximate optically thick radiation 
losses by taking account of reduced escape probabilities in the 
damping wings of strong lines; it is possible that these meth- 
ods can be extended in a robust way to future time-steady stel- 
lar wind models. 

3.3. Heat Conduction 

For the radially oriented flux tubes that we consider, the 
conductive energy exchange rate is dominated by the diver- 
gence of a parallel heat flux density, i.e., 

Qc 0a , = -~( qi \A) . (13) 

Because the solar atmosphere undergoes a transition from be- 
ing strongly collisionally coupled (at low heights and high 
densities) to being nearly collisionless (at large heights and 
low densities), it was realized long ago that the classical 
Spitzer-Harm (SH) prescription for thermal conductivity must 
break down somewhere in the corona and solar wind. We thus 
follow Wang (1993), and others, by using a semi-empirical 
bridging law between the SH heat flux gsH and a completely 
collisionless "free-streaming" heat flux qp§, 

t'coll qm + I'exp ?FS 

n = 1 

t'coll + ^ exp 

where ^ exp = (uj p)\dp/dr\ is the local wind expansion rate 
and v m \\ is the electron-electron Coulomb collision frequency, 



(14) 



In A, 



275 s V10 6 cm- 



10 6 K 



-3/2 



(15) 



lnA ee = 23.2 + -ln [ 

2 \10 6 K 



(Braginskii 1965; Olsen & Leer 1996). The electron Coulomb 
logarithm is approximated by 

7 - llnf * -). (16) 
2 V10 6 cm- 3 / V ; 

A more accurate version of equation (TUfl i was derived by 
Cuperman & Dryer (1985). Under a number of simplifying 
approximations, though (such as nearly Maxwellian distribu- 
tions), their expression reduces to something very close to the 
above bridging formula. 

The radial distance where the heat flux undergoes the tran- 
sition from collisional to collisionless can be estimated by lo- 
cating the point at which v exp = ^coii- For the models of low- 
density high-speed solar wind streams presented below, this 
occurs at about r w 10 solar radii (Rq); for the models of 
high-density low-speed solar wind, this occurs at r m 50-80 
Re- 

In the collisionally dominated limit, we assume SH conduc- 
tion, 

dT 



qsH = -K- 



dr 



(17) 



where the thermal conductivity is assumed to be dominated 
by contributions from free electrons and neutral hydrogen, k = 
FK e + KH (e.g., Nowak & Ulmschneider 1977; McClymont & 
Canfield 1983). The effects of protons and heavy ions are 
neglected because they tend to be overwhelmed by the elec- 
tron conductivity in regions of appreciable ionization (see also 
Ulmschneider 1970; Hansteen & Leer 1995). The electron 
conductivity is given by 

7>5/2 

K e = (1.84x lO^ergcnfV'K" 7 / 2 )— — (18) 

In A eP 
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(Spitzer 1962; Braginskii 1965). Partial ionization effects 
(i.e., F and kh) are parameterized using the results of Shmel- 
eva & Syrovatskii (1973) and McClymont & Canfield (1983), 
with 

l+4.49w + 3.37w 2 + 0.59w 3 

F = (19) 

(l + 3.86w+0.94w 2 ) 2 ' v ; 

where w is the ratio of the electron-electron collision time r ee 
to the electron-neutral hydrogen time r c o, and 



TeO 



n e \nKe \ 52705 K 
Similarly, the neutral hydrogen conductivity is given by 

29.6 T 

K H = 



l + (n p /no)(T/7.6x lO 5 ^" 1 / 2 



(20) 



(21) 



The above partial ionization effects are generally unimpor- 
tant in the solar atmosphere, but we include them for some 
attempt at completeness, and with anticipation that the cooler 
and more optically thick chromospheres of other stars may 
depend more sensitively on them. 

In the collisionless limit we use the free-streaming heat 
flux derived for escaping electrons by Hollweg (1974, 1976), 
which was based on empirical constraints from Forslund 
(1970) and Perkins (1973) on the electron velocity distribu- 
tion, 



<7fs = -a c n e uk B T , 



(22) 



where a c is an order-unity correction factor that depends on 
how the wings of the electron velocity distribution depart 
from a Maxwellian shape. We take a constant value of a c = 4 
as has been often used in solar wind modeling (see, e.g., Leer 
et al. 1982; Scudder & Olbert 1983; Withbroe 1988; Canullo 
et al. 1996; Landi & Pantellini 2003). Tests using Hollweg's 
(1974) more detailed prescription for how a c should depend 
on the solar wind speed and electron temperature produced 
no more than a 10% change in T(r) from the a c = 4 model 
(and no more than a 2% change in the mass loss rate or termi- 
nal wind speed). The simpler constant value was used in all 
subsequent ZEPHYR models. 

Note that equations (fT~4-b and (l22l are in some ways similar 
to the heat flux relations used by Smith & Auer (1980) and 
others for flare plasmas, where the classical heat flux is as- 
sumed to saturate at a value no larger than a threshold flux 
proportional to v\ (where v e is the electron thermal speed). 
The free-streaming heat flux given above is dimensionally 
similar to the saturated heat flux, but the former is propor- 
tional to mv 2 rather than v\. These two limits are related, 
though, because saturation occurs when the electron-electron 
collisional mean free path becomes larger than the local tem- 
perature scale height. When heat can no longer be carried dif- 
fusively by classical conduction, it may be advected directly 
at some characteristic velocity v c , and the heat flux then be- 
comes proportional to v e v 2 . In the low-density supersonic so- 
lar wind, the characteristic speed for, e.g., the total enthalpy 
flux is v c ~ u. To the extent that a c is an order-unity correc- 
tion factor (that depends weakly on v e ), this applies also to the 
heat flux carried in the non-Maxwellian tail of the electron 
distribution. For a high-density laboratory plasma, though, 
Mannheimer & Klein (1975) showed that v c scales directly 
with v e (see also Smith & Lilliequist 1979; Craig & Davys 
1984). 



4. ACOUSTIC WAVES AND SHOCKS 

We include the time-averaged effects of acoustic waves that 
propagate parallel to the magnetic field. The only source of 
these waves that we consider is the convective motion below 
the photosphere which channels compressive wave energy 
into the magnetic flux tubes. Deep in the atmosphere these 
fluctuations take the form of longitudinal, or sausage-mode 
tube waves (e.g., Spruit 1982; Roberts 2000), but because the 
individual flux tubes appear to merge together somewhere in 
the low chromosphere into a region filled with magnetic field 
(Cranmer & van Ballegooijen 2005) we treat these waves as 
standard acoustic oscillations and do not consider thin-tube 
dispersive effects. A separate source of compressive waves, 
which we do not model, may be the gradual parametric decay 
of nonlinear MHD waves in the outflowing stellar wind (see 
Sagdeev & Galeev 1969; Goldstein 1978; Jayanti & Hollweg 
1993). It is usually assumed that these waves have very low 
frequencies and thus are likely to have small rates of damping 
and steepening, and thus a minimal impact on chromospheric 
or coronal heating. 3 

Below we describe how ZEPHYR models the propaga- 
tion, steepening, and dissipation of individual monochromatic 
"packets" of acoustic wave energy (§ 4.1) and how the com- 
plete power spectrum is specified (§ 4.2). 

4.1. Monochromatic Wave Train Evolution 

An arbitrarily steepened acoustic wave/shock train trav- 
els along the field with constant frequency u> and a radially 
varying wavenumber k\\ determined by the dispersion relation 
u> = (u+c s )k\\. The sound speed is given by c 2 s =^fP/p. The en- 
ergy density Us of linear acoustic fluctuations obeys an equa- 
tion of wave action conservation, 



d_ 

dt 



(Us\l_d_ 


~(u + c s )AU s ~ 


\uj'J Adr 


CO* 



(23) 



(e.g., Jacques 1977; Koninx 1992). The Doppler shifted fre- 
quency in the frame of the accelerating wind is given by 
uo 1 = ui-uk\\ and the wave energy density is given by 



1 , 
U s = -pv 
s 



(24) 



where s is a dimensionless shape factor determined by the spa- 
tial profile of the waves (see below). For a small-amplitude 
sinusoid, s = 2, and for a fully steepened sawtooth or N-wave, 
5 = 3. The parallel velocity variance, or squared wave ampli- 
tude, is specified as Vm . Below we give the wave energy flux 
Fs as a lower boundary condition for the ZEPHYR code; this 
quantity is converted into wave energy density using 



(7 + 3) M 



+ c s 



Us 



(25) 



Note that in the limit of a static plane-parallel atmosphere 
(u = and A = constant) equations (f23l>— (t25T> simplify into 
the standard flux conservation quantities implicit in wave- 
heated models of the solar atmosphere since the initial work 
of Schwarzschild (1948), Biermann (1948), and Osterbrock 
(1961). 

3 See, though, Suzuki & Inutsuka (2005) for an example of naturally pro- 
duced compressive waves by similar nonlinear couplings. This kind of pro- 
cess could account for a substantial fraction of the low-frequency density 
fluctuations measured in interplanetary space by in situ spacecraft and radio 
scintillations. 
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The right-hand side of equation (1231 couples the dissipation 
of acoustic wave energy to the Qs heating term in equation (01. 
We include two sources of wave damping: heat conduction 
and entropy gain at shock discontinuities. For the present 
models we ignore radiation damping of the waves (which 
damps high frequencies at low heights, but does not techni- 
cally provide heat to the plasma), viscosity, and ion-neutral 
friction. The acoustic heating rate is given by 



Qs = 27cond£/.S + 



pTAS 

2tt/lo ' 



(26) 



where 7 con d is the linear damping rate due to heat conduc- 
tion and AS is the net entropy jump across a shock (which is 
nonzero only above the height where the wave train has steep- 
ened into shocks). 

The damping rate due to heat conduction is given for adia- 
batic waves (7 = 5 /3) by 

4uj 2 k 

7cond = — -j — (27) 
15kb<7«h 

(e.g., Landau & Lifshitz 1959; Hung & Barnes 1973; Whang 
1997). We use the same thermal conductivity k as is used in 
the classical SH heat flux (eq. fTTIl ). For simplicity, the elec- 
tron and neutral hydrogen conductivities are assumed to dom- 
inate proton and heavy ion heat conduction as well as other 
classical transport processes. Note that 7 con d is the damping 
rate for the wave amplitude; the damping rate for wave energy 
is twice that value (see the factor of 2 in eq. 1261 ). 

The major contribution to the acoustic heating rate comes 
from shock steepening and dissipation. The gain in internal 
energy across an ideal inviscid shock is given by 



TAS 



T 2 -Ti {pi/ pi) 



7-1 



(28) 



where subscripts 1 and 2 denote quantities measured on the 
upstream (supersonic) and downstream (subsonic) sides of 
the shock (Landau & Lifshitz 1959). The above expression 
does not contain the internal energy component from PAV 
work but only the energy that goes into dissipation. Equa- 
tion d26| i uses the approximation from so-called "weak-shock 
theory" that the volumetric heating rate is given by the inter- 
nal energy dissipated at one shock divided by the mean time 
between shock passages in a periodic train. This assump- 
tion breaks down for very strong shocks in the chromosphere, 
which dissipate their energy in a relatively narrow zone be- 
hind the shock (e.g., Carlsson & Stein 1992, 1997), but the 
models presented below do not develop such strong shocks. 

To evaluate equation d28l ) we used the classical Rankine- 
Hugoniot relations for a monatomic (7 = 5/3) gas. These re- 
lations also are valid for a plasma with a constant ionization 
state across the shock, and Carlsson & Stein (2002) found 
that shock trains in the solar chromosphere often approach 
this nearly steady-state condition. The density and tempera- 
ture jump relations can be written in terms of the upstream 
Mach number Mi , 



72 

7] 



pi = ( 7 +l)Mf 
pi (7-l)M2 + 2 

[2 7 MH 7 -l)][(7-l)M? + 2] 
(7+ \fM\ 



(29) 



(30) 



Note, though, that for shocks of arbitrary strength, equation 
(f28b requires the absolute upstream and downstream temper- 
atures Ti and T 2 to be computed, not just their ratio. (In the 



weak-shock limit, T\ w T2 ~ T, the latter being the "mean" 
model atmospheric temperature.) In many astrophysical mod- 
els of shocks, Ti is often assumed to be the undisturbed equi- 
librium state and T2 is computed from equation ( f30b - How- 
ever, in time-dependent simulations of chromospheric shocks, 
T\ is often seen to fall below the time-averaged mean temper- 
ature T and often also below the radiative equilibrium value 
Trad ~ 4500 K. This is believed to arise from adiabatic expan- 
sion behind the shock. 

To compute Ti and T 2 , we first determine the upstream and 
downstream densities p\ and pi relative to the known back- 
ground model density p. With respect to the propagating 
shock train, the background density can be defined as that 
which occurs when the shock passes through zero velocity 
in the reference frame of the undisturbed atmosphere. For an 
ideal sawtooth-shaped N-wave, this occurs for a given height 
at a time halfway between shock passages. Using this defini- 
tion, the analytic results of Bertschinger & Chevalier (1985) 
can be used to estimate the ratio of the minimum (preshock) 
density to the background value to be 



fh _ 1 + Wpi) 
P 2(p 2 /pi) 



(31) 



which then allows both p\ and pi to be computed. To convert 
densities into temperatures, some knowledge of the thermo- 
dynamic cycle of the shock must be incorporated. In other 
words, after the gas is heated and compressed, we need to 
know what "path" it takes as it cools and expands back to 
the preshock values, in order to be heated and compressed 
again as the next shock in the train goes by. Nearly all time- 
dependent models of periodic shocks in stellar atmospheres 
have found that shocks first undergo rapid radiative cool- 
ing at a roughly constant density, followed by nearly adia- 
batic expansion back to the preshock density and temperature 
(e.g., Weymann 1960; Osterbrock 1961; Ulmschneider et al. 
1978; Bowen 1988). This second phase seems to dominate 
the time between shock passages — usually encompassing the 
halfway point used above as the definition for the undisturbed 
density — so we assume that the cooling from the mean tem- 
perature T to the preshock, or upstream temperature T\ is adi- 
abatic, and 

Ti/T = {pi/pV- 1 . (32) 

This, in combination with equations d29l-(l3"Tl>. completes the 
specification of T\ and T2 needed to compute TAS. 

In the limit of low-amplitude shocks (i.e., M\ w 1 + m, 
where m <C 1), equation d28T > reduces to the standard weak- 
shock limit 

27(7-1) , 

(33) 



TAS 



3(7+l> 



3 T 



(e.g., Ulmschneider 1970; Stein & Schwartz 1972, 1973; Mi- 
halas & Mihalas 1984). In the strong-shock limit (Mi > 1 
and pi/pi — > 4), the ratio T\/T approaches a constant value 
of(5/8) 2 / 3 «0.73 andT 2 grows without bound. The first term 
of equation (f28j) dominates the second and 



TAS w c v T 2 « 0.228 c v MJ 



(34) 



The weak and strong limiting expressions are valid to within 
about 25% of the exact result for Mi < 1.2 and M x > 4.6, 
respectively. Because the peak Mach numbers of the shock 
trains in the ZEPHYR models shown below are typically out- 
side these ranges (i.e., M\ 2) we use the full procedure de- 
scribed by equations d28b-d32b, which is valid for shocks of 
arbitrary strength. 
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In order to incorporate the shock dissipation into the acous- 
tic wave transport equation, the wave amplitude vii must be 
converted, where appropriate, into the Mach number M\ . We 
must take account of how the initially sinusoidal wave pro- 
file steepens into a sawtooth shock train. Often shock heating 
is applied only above an estimated shock formation height. 
However, there exists a finite range of heights between the 
first formation of the shock (i.e., where the acoustic wave 
train first "breaks") and the height where the shock train 
has evolved into a complete sawtooth shape. Between these 
heights the velocity amplitude of the shock may only be a 
fraction of the crest-to-trough velocity amplitude of the wave. 

To solve for the wave evolution and dissipation in this re- 
gion, we have computed numerical profiles of both the shape 
factor s and a steepening efficiency factor e. These quanti- 
ties are computed as a function of a dimensionless steepening 
parameter (, which measures by how much a wave crest is 
approaching (or has overtaken) the zero-velocity node imme- 
diately ahead of it. 4 The efficiency e is defined as the ratio of 
the shock velocity amplitude to the full velocity amplitude of 
the (arbitrarily steepened) wave profile, and we define 



1 + 



7+l\ ev 



(35) 



Where the wave is still a sinusoid, e = and Mi = 1, which 
results in no heating. As the wave steepens, e grows from 
to 1. At each radial grid zone, we compute the crest-to- 
node distance factor £ = Az/A, where A is the local parallel 
wavelength and 



Az= T" 



7+1 f v\\dz 



(36) 



where the integration is taken from the lower boundary (at 
which A = Ao) up to an arbitrary height. By computing this 
quantity point-by-point along the grid, the atmospheric strati- 
fication is taken into account accurately. At the lower bound- 
ary, the wave profile is assumed to be a perfect sinusoid, and 
£ = 1/4. This quantity decreases steadily as the wave train 
propagates upward and steepens. When £ reaches zero, the 
shock amplitude has grown to be equal to the wave amplitude 
and the profile is assumed to remain a sawtooth as it propa- 
gates upward. We continue to follow the ever-decreasing £ to 
values below zero, though, because the profile only reaches 
the exact sawtooth "N-wave" shape in the limit of ( — > — oo. 

Figure 2 shows the normalized phase function $(</>) of a 
gradually steepening acoustic wave, where the phase <fi is de- 
fined at a given (constant) height as TT-ujt; it varies from — tt 
to 7r. Figure 2 also plots the efficiency e and the shape factor 
s as a function of £. These have been derived from numerical 
simulations of this steepening sinusoidal wave profile, and we 
give parameterized fits (used by the ZEPHYR code) below. 
The fit to the efficiency is given by 



1 



e(C) = { VWC/Co) 



3.1 







C<o 

o < C < Co 
C>Co 



(37) 



The constant (o denotes the critical breakpoint at which the 
wave train first steepens to infinite slope at the node ahead 
of the crest; it is defined as Co = (tt-2)/4tt 0.091. Note 
that as £ decreases from 0.25 to £o, the efficiency e remains 

4 The node propagates exactly at the linear phase speed, and the crest prop- 
agates faster by a nonlinear factor proportional to the wave amplitude vii . 




phase / 7T 




FIG. 2. — Properties of steepened acoustic waveforms, (a) Normalized 
phase function <3? as a function of the dimensionless wave phase ij> (in units 
of 7r) for a series of steepening parameters f (see labels for values), (b) 
Steepening efficiency ratio e (open diamonds) and scaled shape factor s — 2 
(filled circles) computed numerically for a series of values of Solid lines 
show the analytic fitting functions given in the text. 

zero because no shock transition has yet formed. Only when 
£ decreases below (b does there exist a shock with a finite 
strength. 

The shape factor s is used to convert between wave ampli- 
tude and energy density (eq. l24l ). and it is defined as the 
inverse of the phase-averaged square of the normalized phase 
function, i.e., 



- = ^-l </o[<T>(,.)|- 

s 2n 



(38) 



(see also Koninx 1992; Suzuki 2004). As an example, for a 
sinusoid the average of sin 2 <j> over a full period is 1/2, and thus 
5 = 2. We found the following fit from a series of simulated 
wave profiles undergoing gradual steepening, 



*(0 = 3- 



l + 1.32e" 37 ^ 
l+2.89e-4 4 ^ 



(39) 



(see Figure 2 for a comparison between the numerically de- 
termined values and the fit). For most of the pre-break steep- 
ening (i.e., from £ = 0.25 down to Co) s remains close to 2. 
Note, though, that when the shock grows to full amplitude at 
£ = the shape factor has increased only to about 2.4. The 
additional increase up to 3 occurs for negative values of £. 

4.2. Acoustic Power Spectrum 

The convection zone generates a continuous spectrum of 
acoustic power, and we need to specify this distribution of 
energy as a function of frequency at the photospheric base 
of the model. The power spectrum at larger heights is de- 
termined implicitly as a result of solving the monochromatic 
wave action conservation equations (eq. l23l ) for a range of 
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frequencies. The continuous power spectrum Ps(uj) is defined 
at the base in units of the photospheric acoustic flux, with 



F s = / duP s (u) . 



(40) 



For waves that escape into the upper atmosphere, there should 
be negligible power at frequencies below the acoustic cutoff, 



2c 5 



(41) 



because waves having ui < w ac are evanescent in an ideal 
hydrostatic atmosphere (e.g., Mihalas & Mihalas 1984; see, 
however, Wang et al. 1995; Schmitz & Fleck 1998). We used 
a constant value of uj ac = 0.03 rad s" 1 (i.e., / s» 4.8 mHz, or 
a period of 3.5 minutes), which corresponds to temperature- 
minimum conditions in the upper photosphere. From a theo- 
retical standpoint, this value is slightly on the high side, since 
the photospheric T e s of 5800 K gives a cutoff frequency about 
13% lower. Our higher value for uj dc is an attempt to ensure 
that the modeled frequencies are those that should be prop- 
agating everywhere in the model — i.e., since we do not treat 
evanescent energy loss, we want to include only waves that 
can make it through both the photosphere and the temperature 
minimum region without becoming evanescent. Note, though, 
that from the standpoint of observational helioseismology, our 
value is slightly on the low side (e.g., Jimenez 2006), so it 
seems to be a satisfactory median value. 

The shape of the power spectrum above the cutoff is a sub- 
ject of some ongoing controversy, which we do not attempt 
to address fully. The presence of substantial power at high 
frequencies (f > 20 mHz) is predicted by traditional theories 
of sound generation from turbulent convection (e.g., Lighthill 
1952; Stein 1967; Musielak et al. 1994; Ulmschneider et al. 
1996) and also by observational inferences of time-steady 
chromospheric heating (Kalkofen et al. 1999; Ulmschneider 
et al. 2005; Cuntz et al. 2007). However, evidence also exists 
that there may be an extremely steep decline in the acous- 
tic power spectrum before frequencies of order 20 mHz are 
reached, and thus that high frequencies would not be impor- 
tant to atmospheric heating (e.g., Judge et al. 2003; Fossum 
& Carlsson 2005, 2006). Recent advances in detecting high- 
frequency acoustic fluctuations have been made by Wunnen- 
berg et al. (2002), DeForest (2004), Muglach (2006), and van 
Noort & Rouppe van der Voort (2006), but no firm conclu- 
sions yet exist regarding their impact on chromospheric heat- 
ing. Future observations with higher spatial and temporal res- 
olution are definitely needed. 

Provisionally, we model Ps(w) with a high-frequency tail 
reminiscent of the turbulent convection theories cited above. 
The following parameterization 



Ps(.v) oc 



(w/w max )^/[l + (w/w max ) 2 ^] UJ > UJ d 

UJ < UJ d 



(42) 



has finite power at the cutoff, a peak value at tu mdx > uj ac , and 
a declining tail with P$ oc lo~^ at high frequencies. Typically, 
ip sa 3 and w max is a factor of 2 to 5 larger than w ac . The 
normalization is given by specifying a known value of Fs and 
using equation (l40l . 

Figure 3 illustrates the shape of the acoustic power spec- 
trum that is used in all of the solar models described in this 
paper. We use constant values of w max /w ac = 3 and ip = 2.5. 
The adopted value of ip is the most sensitive to the contro- 
versy over high-frequency acoustic waves, and it deserves 
further discussion. The standard Lagrangian treatment of the 
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FIG. 3. — Normalized power spectra for Alfven waves Pa (solid line) and 
acoustic waves P$ (dashed line). The cyclic frequency / = uj/2tt is given in 
units of mHz, and the period p = l/f is given in units of minutes. 



Lighthill-Stein sound generation mechanism tends to give val- 
ues of ip between 3 and 3.5 in the high-frequency limit for 
both acoustic and longitudinal flux-tube waves (Ulmschnei- 
der et al. 1996; Musielak et al. 2000). However, an alternate 
Eulerian treatment of the turbulent correlations (Rubinstein & 
Zhou 2002) yielded a much shallower decline with increasing 
frequency; ip w 1.3-2.4. We also note that the value ip = 3 
is a special case that corresponds to frequency-independent 
shock steepening for the discrete bins that we use (see below); 
i.e., when ip = 3 the increase in frequency from one bin to the 
next is balanced exactly by the power decrease, such that the 
wave train in each bin steepens into a shock train at the same 
height. We find that a realistic chromospheric temperature 
rise occurs only when higher frequency waves steepen at suc- 
cessively lower heights (which requires ip < 3). Our value of 
2.5 satisfies this requirement while being only slightly lower 
than the range predicted by traditional Lighthill-Stein theory. 

In the ZEPHYR code, the continuous spectrum Ps(ui) is 
modeled as a series of discrete frequency bins, each of which 
is treated independently as described in § 4. 1 . We define the 
bins as "octaves" in frequency space; i.e., the first bin encom- 
passes uj- dc to 2oj. dc , the second encompasses 2ui dc to 4w ac , and 
so on. Ideally, it would have been preferable to use even nar- 
rower bins in order to more accurately represent the shape of 
the spectrum. However, if the frequency bins were too nar- 
row, each would contain a vanishingly small amount of wave 
energy. The calculation of nonlinear steepening (eq. [36)) de- 
pends on the total amplitude that remains reasonably coher- 
ent at a given frequency. Narrow bins — treated independently 
as described above — would essentially destroy this coherence 
and produce an unphysical delay in the onset of steepening 
due to the low amplitudes in each bin. The use of octaves is 
an attempt to balance the needs of frequency resolution and 
realistic coherence for steepening. 5 

5 Too much coherence would also be undesirable. One-dimensional simu- 
lations of chromospheric shocks often result in "shock cannibalization" (i.e., 
overtaking and merging) and an effective filtering out of high frequencies. A 
realistic multidimensional distribution of strong and weak acoustic sources at 
and below the photosphere, though, is more likely to result in an incoherent 
and randomized power spectrum such as we assume here (see, e.g., Cadavid 
et al. 2003; Ulmschneider et al. 2005). 
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In the models computed for this paper, a series of 14 bins 
is used. The power missed by ignoring frequencies above the 
maximum value of 2 14 w ac is only ~ 10~ 5 of the total flux. The 
monochromatic frequencies used in the wave action conser- 
vation equation are computed as mean, or first-moment fre- 
quencies over the power spectrum of each bin. The first one 
(w ac to 2w ac ), for example, has a mean of 1 .62w ac (i.e., / = 7.7 
mHz). 

5. ALFVEN WAVES AND TURBULENCE 

We model transverse MHD fluctuations in the radially 
oriented magnetic flux tube as ideal incompressible Alfven 
waves. The equations and assumptions we use generally fol- 
low Cranmer & van Ballegooijen (2005), who computed var- 
ious properties of Alfven waves from the photosphere to the 
interplanetary medium assuming a known background plasma 
state (m, p, Bq). Here we compute the wave properties and 
background plasma properties together in a self-consistent 
manner. The Alfven waves are believed to arise from trans- 
verse jostling of magnetic flux tubes in the intergranular pho- 
tosphere by convective motions on granular (^1000 km) hor- 
izontal scales. Note that, similar to the acoustic waves, the 
Alfven waves are expected to have modified dispersive prop- 
erties in the photosphere and low chromosphere (e.g., they 
become kink-mode tube waves in regions where the flux tubes 
are surrounded by field-free regions; see Spruit 1981, 1982), 
but we ignore these effects because they occur over a small 
range of heights. Our treatment of Alfven wave reflection and 
turbulent dissipation is complementary to the recent work of 
Verdini et al. (2005, 2006). 

As in the previous section, we first discuss the overall 
wave energy transport and dissipation mechanisms (§5.1) and 
then describe the power spectrum of the fluctuations (§ 5.2). 
The way that the frequency dependence is assimilated into 
the model is slightly different from the acoustic-wave case, 
though. 

5.1. Alfven Wave Action Conservation 

For transverse incompressible waves, we solve a single 
wave action conservation equation with a damping term that 
contains information about the power spectrum and wave re- 
flection. We believe that a monochromatic treatment of the 
wave action conservation (i.e., to treat each frequency bin 
independently from the rest) is inappropriate to use when 
the dominant dissipation mechanism is MHD turbulence. A 
turbulent cascade inherently contains strong nonlocal "mix- 
ing" in frequency and wavenumber space. Thus the turbulent 
damping rate is computed only using quantities that have been 
integrated over the power spectrum (see below). 

We assume the Alfven waves are dispersionless with a 
phase speed given by irf/k\\ = u + Va in the stationary refer- 
ence frame of the Sun. The Alfven speed is defined as Va = 
B /(4tt/9) 1 / 2 . Although we solve the full non-WKB transport 
equation to obtain the relative distribution of outward and in- 
ward propagating waves, we use the following simplified ver- 
sion of the wave action conservation equation to compute the 
radial evolution of the total (frequency-integrated) wave en- 
ergy density Ua- 



d_ 

di 



Ua 



]_d_ 

/ A dr 



(u + V a )AUa 



Qa 



^ (43) 



(see, e.g., Jacques 1977; Isenberg & Hollweg 1982; Tu & 
Marsch 1995). The presence of the Doppler shifted frequency 



u/ (measured in the solar wind frame) in the above equation 
is deceptive, since the wave action conservation is essentially 
frequency-independent. Note that u>' = to—uk^ = u>Va/(u+Va), 
and thus a constant factor of w can be factored out of each 
term in equation d43l . The total Alfvenic energy density Ua is 
the sum of the energy densities in kinetic and magnetic fluc- 
tuations, 



U K = 



U B = 



8vr 



(44) 



where vj_ and B± are the transverse velocity and magnetic 
field oscillation amplitudes. We assume equipartition between 
the two components (Uk = Ub), such that Ua = 2Uk = pv\. 

The use of equation d43b contains the assumption that the 
energy balance is dominated by outward-propagating waves, 
whose energy density U- exceeds the energy density of 
inward-propagating waves U+. By definition, Ua = U- + U+, 
and 

pZl 

V± = L f- (45) 

where the Elsasser (1950) amplitudes are defined here as 
Z± = |vj_ ±B^/(47T/9) 1 / 2 |. For future reference, we give the 
expression for the net outward Alfven wave flux, 



u(U k + 2Ub) + V a (U--U + ) 



(46) 



(e.g., Heinemann & Olbert 1980). The two assumptions of 
kinetic/magnetic equipartition and outward-wave dominance 
in equation d43l tend to break down in the lower atmosphere 
when there is strong non-WKB wave reflection, but in the 
corona — where the dominant Alfvenic heating occurs — these 
assumptions seem to be appropriate (see also Cranmer & van 
Ballegooijen 2005). 

The only physical source of Alfven wave damping that we 
include is the turbulent cascade, which ultimately must termi- 
nate in an irreversible conversion of wave energy into heat. 
This is, in some sense, a controlled experiment to evaluate to 
what degree turbulence may be the dominant cause of coro- 
nal heating and solar wind acceleration, but other damping 
mechanisms for Alfven waves have been suggested. Colli- 
sional damping mechanisms for MHD waves have been stud- 
ied extensively in the context of the dense coronal base (e.g., 
Alfven 1947; Kuperus et al. 1981; Narain & Ulmschneider 
1990, 1996; Porter et al. 1994; Roberts 2000), but the open 
field lines that feed the solar wind have lower densities and are 
thus less collisionally dominated than closed loops in the low 
corona. Cranmer & van Ballegooijen (2005) investigated lin- 
ear viscous dissipation of Alfven waves along an open mag- 
netic flux tube and found that viscosity should be negligible 
for waves having periods longer than about 1 minute. 

We adopted the following phenomenological form for the 
MHD turbulence damping rate, 



Qa = p£\ 



turb " 



ZiZ + + Z\Z- 



4L 



(47) 



(Hossain et al. 1995; Zhou & Matthaeus 1990; Matthaeus et 
al. 1999; Dmitruk et al. 2001, 2002). The transverse length 
scale L±(r) represents an effective perpendicular correlation 
length of the turbulence for the largest "driving" eddies. We 
used the standard assumption that L± scales with the trans- 
verse width of the open flux tube; i.e., that it remains propor- 

-1 12 

tional to B (see also Hollweg 1986). Ideally, the evolution 
of L± should be coupled to the radial variation of the fluctu- 
ation energy (see, e.g., eq. [3] of Matthaeus et al. 1999) as 



CORONAL HEATING AND SOLAR WIND ACCELERATION 



11 



well as to the non-varying field strength. Future work should 
include both this effect and additional tests of whether the 
above phenomenological form is an adequate representation 
of the true anisotropic cascade. 

The normalization of L± is one of the few free parame- 
ters of our model. We note that Cranmer & van Ballegooijen 
(2005) found that L± should be about 1 100 km at the height in 
the low chromosphere where thin flux tubes merge with one 
another. Mapping this value down to the photosphere would 
yield a lower boundary value L± Q = 320 km. This spatial 
scale is intermediate between the probable horizontal size an 
individual flux tube in the photosphere (50-100 km) and the 
size of a convective granule (1000 km). This value is also 
similar to the width of an intergranular lane and also the mean 
separation between photospheric flux tubes in the quiet-Sun 
supergranular network (350-700 km; see Cranmer & van Bal- 
legooijen 2005). 

The factor £ t urb in equation d47l i is a turbulent efficiency 
that accounts for regions where the turbulence may not have 
time to develop before the waves or the wind carry away the 
energy (see Dmitruk & Matthaeus 2003). We estimated this 
efficiency factor to be 



c-turb 



1 



l+(feddy/fref)' ! 



(48) 



where the two time scales in this expression are f e ddy, a non- 
linear outer-scale eddy cascade time, and f le f, a timescale for 
macroscopic Alfven wave reflection. In most of the models 
presented below we take n = 1, based on analytic and numeri- 
cal models of Dobrowolny et al. (1980), Matthaeus & Zhou 
(1989), and Oughton et al. (2006). Dmitruk & Matthaeus 
(2003) found that the turbulent cascade has sufficient time to 
develop and heat the plasma only when f e dd y <C f re f • Thus, our 
efficiency factor above quenches the turbulent heating when 
feddy S> fref, i.e., when the Alfven waves want to propagate 
away much faster than the cascade can proceed "locally." The 
reflection time is defined simply as f le f = 1/| V • Va|, and the 
eddy cascade time is given by 



feddy — 



(l+M A )v A 



(49) 



where the Alfvenic Mach number Ma = u /Va and the numer- 
ical factor of \/3tt comes from the normalization of an as- 
sumed shape of the turbulence spectrum (see Appendix C of 
Cranmer & van Ballegooijen 2005; see also Higdon 1984; 
Shebalin et al. 1983; Goldreich & Sridhar 1995, 1997; Bhat- 
tacharjee & Ng 2001; Cho et al. 2002). When n = 1, the ef- 
ficiency factor provides an approximate bridging between a 
Kolmogorov (1941) scaling, when f e dd y <C f re f, and an IK-like 
(Iroshnikov 1963; Kraichnan 1965) scaling, when feddy 3> fref- 
There is still some controversy, though, over which type of 
cascade rate is appropriate for MHD turbulence in the solar 
wind. 

We separated the Alfven wave energy into outward (Z_) 
and inward (Z+) components by solving a modified form of 
the non-WKB transport equations of Heinemann & Olbert 
(1980), Barkhudarov (1991), Velli (1993), and Orlando et 
al. (1996). 6 These frequency-dependent equations are dis- 

6 By WKB (Wentzel, Kramers, Brillouin) we do not refer to any specific 
asymptotic expansion. The WKB limit of pure outward propagation, with no 
reflection, is amenable to the standard "eikonal" approximation by defining 
a local wavenumber. The treatment of non-WKB reflection is more general 
in that the radial parts of the Z± eigenfunctions are computed numerically 
without the use of a wavenumber. 



cussed below in § 5.2 and their solution provides a spectrum- 
averaged value of the effective local reflection coefficient 
1Z = Z+/Z-. Knowing 1Z and Ua (from eq. 1431 ) allows the 
Elsasser amplitudes to be computed at all heights, 



Z_ = 



4U A 



pd+n 2 ) 



z + = n\z.\ , 



(50) 



and the turbulent damping rate in equation d47l ) is then speci- 
fied. 

5.2. Alfven Wave Frequency Dependence 

For a series of Alfven wave frequencies u>, we solved 
the non-WKB wave transport equations in the dimensionless 
form given by equations (24) and (33) of Barkhudarov (1991), 

_ (* 2 -i)cosr 

dr 



2H A 

dT _ 0J 2 +l)sinr_ 2uV A 
~dr~ ~ 



(51) 



2H A V 



V 



(52) 



where ^ is related to the frequency-dependent reflection co- 
efficient IZu via 



u-Va 
u + Va 



(53) 



T is the angular phase shift between the Z_ and Z+ wave trains, 
and Ha is the (signed) scale height for the Alfven speed, or 
V A /(dV A /dr). Although the models of Barkhudarov (1991) 
were limited to spherical geometry (A cx r 2 ) the above rela- 
tions are valid for any A(r) (see also Cranmer & van Bal- 
legooijen 2005). We followed the general solution proce- 
dure outlined by Barkhudarov (1991) for integrating across 
the Alfvenic singular point ta, where u = Va and thus '3/ = 0. 
The reflection coefficient remains finite at this point, and it 
can be written exactly as 



sjuj 2 + (u'- 



VZ) 2 



(54) 



where u' and V A are radial derivatives of the outflow speed 
and Alfven speed taken at r = rA ■ The ZEPHYR code utilizes 
fourth-order Runge-Kutta numerical integration to solve for 
Huj(r) above and below r A . We ensured that \TZ U \ < 1 at all 
heights and for all frequencies. 

The dissipation of Alfven waves is explicitly not included 
in the equations given above for the reflection coefficient. Al- 
though it is possible to include nonlinear damping consistent 
with equation d4Tb in these transport equations (see, e.g., Ver- 
dini et al. 2005, 2006), we remain cautious about the combina- 
tion of strong turbulent damping and "monochromatic" wave 
quantities. Additional simulations may be required in order 
to better guide the use of phenomenological nonlinear terms 
when following the development of a spectrum of Alfvenic 
fluctuations. 

The full frequency-averaged reflection coefficient is com- 
puted by weighting IZ^ by the Alfven wave power spectrum 
Pa(oj), with 

/ dajP A (uj)TZ 2 Jr) 



J dwP A {to) 



(55) 



where the square of 1Z is used because the power spectrum 
is an energy density quantity and 1Z is a ratio of amplitudes. 
The spectrum is used essentially as a weighting function, and 
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we assume it has a constant shape as a function of height (see, 
e.g., Figure 8 of Cranmer & van Ballegooijen 2005). 

Figure 3 shows the shape of the Alfven wave power spec- 
trum that we use in the ZEPHYR models discussed below. 
Specifically, this is the power spectrum of total (kinetic plus 
magnetic) energy from the model of Cranmer & van Balle- 
gooijen (2005) taken at the height of the transition region. It 
was computed from an empirically constrained photospheric 
spectrum of the dynamics of thin flux-tubes, which is de- 
scribed as a linear superposition of two types of motion. First, 
isolated flux tubes undertake "random walks" in response to 
convective granulation; they have a power spectrum that has 
been constrained by observed G-band bright point motions 
(e.g., van Ballegooijen et al. 1998; Nisenson et al. 2003). 
Second, flux tubes exhibit sporadic rapid horizontal "jumps" 
that probably represent merging, fragmenting, or reconnect- 
ing with surrounding magnetic fields (e.g., Choudhuri et al. 
1993; Berger & Title 1996; Berger et al. 1998; Hasan et al. 
2000) which we modeled as a series of periodic impulsive 
motions. High-resolution photospheric observations thus pro- 
vided the kinetic energy spectrum, and the partitioning be- 
tween kinetic and magnetic energy components (needed to 
compute the total energy spectrum) was determined using the 
analytic linear theory of kink-mode waves in a stratified atmo- 
sphere. Finally, propagation effects between the photosphere 
and the transition region (z = 0.003 Rq in the model of Cran- 
mer & van Ballegooijen 2005) were taken into account by 
solving the combined non-WKB kink-mode and Alfven-mode 
transport equations for an empirically constrained (FAL-C') 
background plasma state. 

There are several noteworthy features of the spectrum Pa{oj) 
in Figure 3. There is power at the lowest frequencies, corre- 
sponding to periods of hours to days, because of the assumed 
random-walk component of the photospheric flux tube mo- 
tion (with Pa oc e~ UT , where r w 60 s). However, the rela- 
tive amount of power at periods longer than 1 hour is negligi- 
ble compared to the total; this stands in contrast with in situ 
measurements that show the majority of power to be at these 
long periods (e.g., Goldstein et al. 1995; Tu & Marsch 1995). 
There is growing evidence that the low-frequency fluctuations 
seen in interplanetary space may be the result of the passage 
of multiple uncorrelated flux tubes past the spacecraft and 
not intrinsic turbulence within any one flux tube (see, e.g., 
McCracken & Ness 1966; Jokipii & Parker 1969; Bruno et 
al. 2001; Giacalone & Jokipii 2004; Giacalone et al. 2006; 
Borovsky 2006). 

The damping of evanescent kink-mode waves can be seen in 
Figure 3 from the mild discontinuity at the cutoff frequency 
of 1.4 mHz. Only about 40% of the energy of evanescent 
waves is lost, though, because the flux tubes remain thin (and 
thus gravitationally buoyant) only below a low-chromosphere 
"merging height" of 600 km. Above this height the magnetic 
fields expand horizontally to fill the volume above supergran- 
ular network lanes and the transverse waves both below and 
above the cutoff can propagate freely as Alfven waves. The 
multiple peaks in the spectrum between about 3 and 20 mHz 
are a result of propagation effects between the photosphere 
and the merging height. In the model of Cranmer & van Bal- 
legooijen (2005) the atmosphere was modeled with the non- 
isothermal FAL-C' temperature structure and a radially vary- 
ing filling factor for the flux tubes; both of these factors re- 
sulted in a complicated frequency dependence for the non- 
WKB transmission of kink-mode waves between z = and 
600 km. Coincidentally, the frequencies of the maxima in 



both Pa and P$ are each about a factor three higher than their 
respective cutoff frequencies. 

It is important to note that the Alfven wave amplitude vj_ 
(obtained from the WKB-like eq. l43l ) diverges from the ac- 
tual transverse velocity of oscillating magnetic flux tubes in 
the lower atmosphere. The effects of evanescence and non- 
WKB wave reflection are not directly included in equation 
d43b . but they end up being of minimal importance in the 
corona and solar wind. We thus can estimate the true velocity 
amplitude w± after the final iterated solar atmosphere param- 
eters have been determined, as 

l+Hj \ { 1 > ^ > ^kc 

\-TZl ) X 1 exp 1 - yj\-{u}/u) kc ) 2 , uj<uj kc ' 

(56) 

where uj kc is the kink-mode cutoff frequency. Above the mid- 
chromosphere height where thin flux tubes merge with one an- 
other, cjkc effectively goes to zero and there is no evanescence. 
The IZu, factor above in parentheses corrects for the approxi- 
mation that equation (l43l follows only outward-going waves. 
For the ZEPHYR models presented below, the spectrum- 
averaged value of w± at the photosphere is typically a factor 
of 2 to 5 times larger than vj_ at this lower boundary. Above 
the transition region, though, w± w vj_ (see also Figure 1 1 of 
Cranmer & van Ballegooijen 2005). Despite the fact that vj_ 
underestimates the actual velocity amplitude, we believe it is 
more appropriate to use this as an imposed lower boundary 
condition rather than to use w± . The latter quantity depends 
on quantities that are computed self-consistently along with 
the other plasma parameters and are not known a priori. 

6. WAVE PRESSURE ACCELERATION 

Just as electromagnetic waves carry momentum and exert 
pressure on matter, propagating acoustic and MHD waves can 
also do work on the mean fluid via a similar kind of radia- 
tion stress. (Bretherton & Garrett 1968; Dewar 1970; Belcher 
1971; Alazraki and Couturier 1971). For parallel-propagating 
acoustic and Alfven waves, the time-averaged radial wave 
pressure acceleration was derived in detail by Jacques (1977) 
to be 

1 dU A fl+l \ dU s _U S OA 
dr 



D = 



(57) 

2p dr \ 2p J dr Ap dr 

As above, we ignore departures from kinetic-magnetic energy 
equipartition for the Alfven waves; this ends up being a good 
approximation for solar wind conditions (see, e.g., Heine- 
mann & Olbert 1980; Cranmer & van Ballegooijen 2005). 

To implement the above expression in the ZEPHYR code, 
we used the wave action conservation equations d23b and d43l 
to reformulate the momentum conservation equation into a 
time-independent Parker (1958, 1963) critical point equation 
with additional terms that affect the definition of the critical 
point. Assuming 7 = 5/3 throughout, the modified momen- 
tum equation becomes 

■ 2X du GM* , , d\\\A 
dr 



, 9 in dlna 2 1 
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dr 
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dr 
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(58) 



2(u + V A ) 3(u + c s )_ 

where a is the isothermal sound speed (k^T /mu) 1 ^ 2 and the 
modified critical speed is given by 

U A 
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with the bulk-flow Mach numbers Ma = u/Va and Ms = u/c s . 
In the presence of acoustic waves, the additional terms above 
are given by 



Us f 1-7 Ms 



3p \ \ +M s 
4Us (Ms-\ 



3p \M s +l 



(60) 



(61) 



The modified critical radius r c is found by locating the point 
where the right-hand side of equation ( l58l is zero, and thus 
u = u c . 

For rapidly expanding superradial flux tubes, Kopp & 
Holzer (1976) pointed out that there may be more than one 
possible location for the critical point. Defining the right- 
hand side of equation (IB"8l as the radial derivative of a known 
function J-(r), they essentially found that the global minimum 
of T(r) specifies the critical point location that allows for a 
stable and continuous solution for u(r) from the Sun to in- 
terplanetary space (see also Vasquez et al. 2003; Cranmer 
2005a). The Kopp & Holzer (1976) result, though, was de- 
rived for an ideally polytropic corona without wave pressure 
acceleration. It is not clear if this global minimum condition 
is applicable to the models described in this paper. For all 
of the ZEPHYR models discussed below in § 8, though, the 
global minimum of T(r) coincides with the largest possible 
value of r c . For expediency, then, when there are multiple 
possibilities we choose the largest value for the critical radius 
that also exhibits a minimum in T (i.e., d 2 J- '/dr 2 > 0). 

Note that equation ( f58l > includes the wave pressure acceler- 
ation, but not the explicit derivatives of Us and Ua that appear 
in equation ( 1571 ). The modified momentum equation thus alle- 
viates the need to perform noisy numerical differentiation (see 
also Jacques 1977; Hartmann & MacGregor 1980; DeCampli 
1981; Holzer et al. 1983; Wang & Sheeley 1991). The damp- 
ing rates Qs and Qa appear explicitly in the momentum equa- 
tion, which provides additional nonlinear coupling between 
the momentum and energy equations. 

Wave pressure acceleration has been invoked by Laming 
(2004) as a potential explanation for the relative enhancement 
of low FIP (first ionization potential) elements relative to high 
FIP elements in the corona and solar wind. 7 If Alfven waves 
exert an appreciable force on ions in the chromosphere and 
transition region (where there are still neutrals that do not feel 
this force) significant "fractionation" may occur. We apply 
Laming's (2004) idea to the models presented below by eval- 
uating his integrated momentum equation for atoms and ions 
of element s undergoing fractionation. A slightly simplified 
version of this equation (which assumes that there is no ex- 
plicit dependence of D on the density gradient) is 



In 



(PsV s )u 
(PsV 2 )l 



dz- 



D 
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(62) 



where p s and v s = (c 2 + v 2 p 1//2 are the element's mass density 
and effective parallel turbulent speed. The ionization fraction 
of the element — essentially one minus the neutral fraction — is 
given by and the collision rates between ions and neutrals 
of element s and the ambient gas are given by v s i and v s n (see 
Laming 2004 for detailed expressions). Subscripts I and u de- 
note quantities computed at heights z\ (a lower boundary we 

7 For other potential explanations of the FIP effect, see also von Steiger 
& Geiss (1989), Vauclair (1996), Arge & Mullan (1998), Schwadron et al. 
(1999), and references therein and in Laming (2004). 



take as the photosphere) and z u (an upper boundary that can 
be anywhere in the corona or solar wind), respectively. The 
degree of FIP enhancement between these heights is obtained 
by solving equation d62l twice: once for the ratio p s , u / p s ,i for 
a low FIP element, and once for a high FIP element. Dividing 
one ratio by the other cancels out the overall density strat- 
ification between zi and z u and leaves only the fractionated 
abundance difference. 

7. SOLUTION METHOD 

We compute time-independent solutions to the conserva- 
tion equations given above by applying a new hybrid method 
of iteration and relaxation. Reasons for not following time- 
dependent fluctuations explicitly were summarized in § 2. It 
has become common, though, to use time-dependent hydro- 
dynamics codes to find stable time-steady solar wind solu- 
tions (see recent work by, e.g., Li et al. 2004; Lionello et al. 
2005; Lie-Svendsen & Esser 2005, and references therein). 
We break from this tradition for several reasons. First, in 
the one-fluid case, the time-independent mass and momen- 
tum conservation equations are relatively easy to solve, and 
only the energy equation requires special treatment. Second, 
when modeling the photosphere, chromosphere, corona, and 
wind all in one grid, the energy equation changes its basic 
character at different locations depending on which terms are 
dominant. In the corona, the Spitzer-Harm conduction makes 
it a second-order parabolic differential equation. In the outer 
solar wind the gradual transition to collisionless heat conduc- 
tion reduces it to a first-order differential equation. In the 
photosphere and chromosphere it is essentially a zeroth-order 
differential equation — i.e., an algebraic balance between heat- 
ing and cooling. These changes, combined with the huge dy- 
namic range in quantities such as density, temperature, and 
wave phase speeds, make it difficult to implement and opti- 
mize a robust time-dependent numerical scheme. Even im- 
plicit methods, which are not necessarily limited by propaga- 
tion across the smallest grid zones, appear to be prohibitively 
difficult to set up properly when the source terms (i.e., the 
waves) have such a complex nonlinear dependence on the pri- 
mary plasma parameters. 

The code developed for this work is called ZEPHYR. Af- 
ter setting up an initial trial guess for the plasma parameters 
p, u, and T, the code iterates a fixed number of times (typi- 
cally 100 to 200) alternately between solutions of the energy 
equation (eq. ifJJ) and the other constitutive and conservation 
equations (eqs. ffl, 0, g), 0, 0, ED, ED, ED, (52)) 
to find a steady-state accelerating wind solution. To stabilize 
the iteration process and avoid pathological solutions, most of 
the plasma parameters are "undercorrected" using a scheme 
described below — i.e., rather than replacing the old solution 
with the new one, a fractional step toward the new solution is 
taken. 

The spatial grid used by ZEPHYR has variable zone spac- 
ing that depends on the height above the lower boundary. The 
total number of grid zones N is divided into two subsets: 45% 
of the zones are allocated to a fine mesh with constant spac- 
ing in the lower atmosphere (i.e., between z = and a fixed 
zone-midpoint height z m id = 0.005 Rq) and 55% of the zones 
occupy the rest of the grid with spacing that increases by 1 % 
per zone from the value at z m id up to the top of the grid. All 
models described in this paper have N = 1300 and a fixed grid 
spacing of 8.56 x 10~ 6 ^ Q w 6 km in the lower atmosphere. 
At the top of the grid (z = 1200i? Q » 5.6 AU) the grid spac- 
ing has increased to 24 R & . The relative spacing Ar/R Q thus 



14 



CRANMER, VAN BALLEGOOIJEN, & EDGAR 



increases from about 10~ 5 in the lower atmosphere up to 0.02 
at the top of the grid. 

The initial condition for the temperature is a gray radia- 
tive equilibrium atmosphere (i.e., T = T m & as defined in § 3.2) 
that is perturbed by transitioning (using a hyperbolic tangent 
function) to a corona/wind value above a specified transition 
region height of 0.9z m id. The temperature above this transition 
obeys a power law in r, with a basal value of 1 .5 x 10 6 K and 
a very slow radial decline proportional to f~° 05 . The lower 
boundary (z = 0) is defined as tr = 2, which corresponds to 
a base density p© = 2.45 x 10~ 7 g cm" 3 . This choice for the 
lower boundary condition is used in order to include the tra- 
ditional photosphere (either tr = 1 or 2/3) in the interior of 
the grid. The temperature at the lower boundary is defined by 
the gray atmosphere condition (eq. fTOl ). The initial density 
distribution is assumed to be hydrostatic (i.e., u = 0). Refined 
initial guesses for the density and outflow speed are evalu- 
ated as the first steps in the iteration process described below. 
Tests have shown that finding the proper solution is insensi- 
tive to the details of the initial guess, but the iteration method 
converges faster when the initial guess has properties closer 
to a realistic solar atmosphere. 

The main "outer" iteration loop consists of two interior 
modules, each of which undergoes a number of "inner" itera- 
tions. The first module solves for the dynamics (p, u) and for 
the various source terms in the momentum and energy equa- 
tions. The five steps taken in one inner iteration of this module 
are as follows. 

1 . The hydrogen ionization fraction x, the gas pressure P, 
and the internal energy density E are computed as de- 
scribed in § 3. 1 . The optical depth scale tr is integrated 
and the radiative cooling/heating rate Qiad is computed 
(§ 3.2). The heat conduction rate Q con d (§ 3.3) is 
also determined using the four-point finite differenc- 
ing scheme discussed below for the radial derivatives 
dT/dr and dqn/dr. The heat conduction Qcond is ar- 
tificially suppressed in the outermost 10 grid zones in 
order to produce a more well-behaved upper boundary 
condition that is dominated by the outward advection of 
all characteristics. (These 10 zones are not considered 
to be part of the actual solar wind solution.) 

2. The modified Parker critical point equation d58l l is 
solved for u(r) by integrating up and down from r c . To 
step from the critical point (which generally falls be- 
tween grid zones) to the grid zones immediately above 
and below, an analytic derivative (du/dr) c computed 
from L'Hopital's rule is used in order to avoid the well- 
known instability at an X-type singular point. At all 
other grid zones, fourth-order Runge-Kutta integration 
is used (e.g., Press et al. 1992). The undercorrection 
scheme described below is used for u(r). 

3. The time-steady mass conservation equation (Q]i is 
solved for p(r) in a straightforward way. The base den- 
sity pQ described above is kept fixed and the prior step's 
solution for the basal outflow speed m q = u(R Q ) is used 
to determine the constant mass flux puA. The density 
is then computed exactly at each grid zone, but the un- 
dercorrection scheme is also used for p(r) in order to 
prevent too rapid a change. 

4. The properties of the acoustic wave spectrum are com- 
puted at each grid point and the spectrum-integrated 



values of Us and Q$ are determined (§4). The steep- 
ening of each monochromatic wave train is computed 
by integrating equation (f36t up from the lower bound- 
ary simultaneously with the wave action equation ( l23l . 
Simple first-order Euler steps are used to perform both 
integrations. 

5. The non-WKB Alfven wave transport equations (§ 5.2) 
are solved for a number of frequency dependent re- 
flection coefficients 1Z u {r) which are summed over 
the power spectrum to obtain IZ(r). The adaptive 
fourth-order Runge-Kutta scheme developed by Cran- 
mer & van Ballegooijen (2005) has been included in the 
ZEPHYR code. Once the Z± Elsasser variables have 
been determined, the heating rate Qa is computed and 
the Alfven wave action conservation equation is inte- 
grated up from the lower boundary to obtain £/a (§ 5.1). 

These steps are repeated for a fixed number of inner iterations 
(typically 50 to 100) in order to reach internal consistency. 
However, because the non-WKB wave reflection equations 
( 15114521 dominate the computation time, the reflection coeffi- 
cient 1Z is recomputed only for the first two of these iterations. 

The second main module of ZEPHYR solves the energy 
conservation equation and obtains a time-independent tem- 
perature distribution T(r). There are three main steps that are 
repeated until either a certain degree of convergence has been 
attained ((6E) < 10" 3 ) or a maximum number of iterations 
(typically 1000) is exceeded. The convergence quantity (5E) 
is defined to be an average over the relative convergence at 
each grid zone. For the radial grid with N discrete zones, 

W^i-YfJ&L-) (63) 

where i denotes quantities computed at each grid zone. The 
numerator is obtained by solving equation (0 for the explicit 
time dependent term E = dE jdt at each grid zone and taking 
the absolute magnitude. This residual-like quantity is nonzero 
for solutions that have not yet converged to a steady state. The 
denominator is the absolute magnitude of the largest single 
term in equation ©; this includes the two advection terms on 
the left-hand side, and it also includes the separation of the 
optically thick J and S terms in Qthick (eq- 0) which balance 
one another in the photosphere. The solution that satisfies ex- 
act time-independence (for the energy equation) would have 
(SE) = 0. The following three steps repeated by the second 
module of ZEPHYR are designed to hone in on this solution. 

1. The core procedure in solving the energy equation is 
relaxation using E at each grid zone. The sign of E 
is used to determine whether the current solution for T 
should be increased or decreased, and the magnitude of 
the change is computed from a positive-definite correc- 
tion factor c(r). 8 This factor is initialized to a constant 
value of 0. 16 at the start of the inner iteration loop. Dur- 
ing each relaxation step it is either kept constant, if E 
has kept the same sign from the last step to the cur- 
rent step, or it is reduced in magnitude by 7%, if £ has 
switched signs from the last step to the current step. The 
reduction in c that occurs when E oscillates in sign is 

8 The fact that the magnitude of E is not used is the main reason that this 
technique is called "relaxation" and is not really a variety of time-dependent 
hydrodynamic evolution. 
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a kind of "annealing" that allows the relaxation method 
to focus in on the time-independent solution (i.e., the 
solution with E = everywhere). The updated tempera- 
ture at each grid zone is thus determined by multiplying 
the old value by a factor of (1 +cE/\E\). The optimized 
numerical values given above were determined from a 
large number of tests with a range of values. 

2. Because the relaxation method above can lead to dis- 
continuous jumps in T(r), we perform trial piecewise 
smoothing in order to reduce unphysical fluctuations. 
The radial grid with N = 1300 is broken up into 130 
pieces each having 10 zones. (At the start of each itera- 
tion, the offset point that defines the start of the first 10- 
zone piece is shifted by one grid zone.) If a piece con- 
tains more than one change of sign in the slope dT/dr, 
it is smoothed using a Gaussian filter (with a two zone 
half-width) until only or 1 changes of sign remain. 
The piecewise nature of this smoothing is necessary so 
that the entire grid does not "suffer" when just a small 
part of it contains numerical noise. 

3. The convergence parameter (SE) for the current inner 
iteration is compared to the best solution (i.e., the low- 
est value of (SE)) that has been found during this outer 
iteration loop. If the solution has improved, the best 
solution is updated. If the solution has gotten worse, 
we discard it and revert to the saved best case. Note, 
though, that this comparison is not performed at every 
inner iteration step. (Doing it every time could lead to 
an infinite loop with no changes ever made to T.) Tests 
showed that it is best to allow the solutions to evolve 
for a while and only perform this comparison every 15 
to 20 iteration steps (we use 17). 

The converged value of T that emerges from this module is 
undercorrected before starting the next outer iteration. 

The undercorrection scheme that is used at various points 
in the ZEPHYR code was motivated by globally convergent 
backtracking methods for finding roots of nonlinear equations 
(e.g., Dennis & Schnabel 1983). Rather than taking the full 
suggested iteration step, which may propel the solution away 
from the desired region of convergence, it is sometimes best 
to take only a partial step. For a scalar quantity fjj at radial 
grid zone i and iteration step j, we specify this partial step as 



fij+i - fij 



fij+i 



fu 



(64) 



where fij+i is the next suggested iteration that was obtained 
by solving one of the conservation equations. The exponent e 
describes the degree of undercorrection. When the solutions 
are nearly converged, the full iteration step should be taken 
(f«l). When the solutions are far from convergence, though, 
we require substantial undercorrection (0 < e <C 1). To obtain 
this exponent, we use 



e = e + (l-eo)min 



fj 



fj+i 



(65) 



where the minimum is taken over the entire radial grid (thus 
the absence of ;' subscripts) so that the worst agreement be- 
tween the current iteration and the next suggested iteration 
is highlighted. Because fj +l may be larger or smaller than 
fj, both the ratio and its reciprocal are used above, and the 
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FIG. 4. — Convergence of the ZEPHYR outer iteration process toward a 
steady-state solution. From top to bottom, the plasma parameters shown in- 
clude the solar wind speed Uoo at the top of the grid (dotted line), the spher- 
ical mass loss rate M (upper solid line), the maximum coronal temperature 
(dashed line) — all in scaled units as listed in the captions — and the dimen- 
sionless energy convergence parameter (SE) (lower solid line). 

largest value that the minimum can take is 1. Tests have 
shown eo = 0.17 to be a robust value; it represents a practical 
lower limit to e that alleviates making infinitesimally small 
corrections. 

When radial derivatives need to be taken in the ZEPHYR 
code, we use the following four-point finite differencing 
scheme, 



df 
dr 



= C 



fi + l-fi-l\ +{l _ Q ( fi + 2-fi-2 \ (66) 



r i+ 2-n-2 



for discretized quantities on the radial grid that has un- 
equally spaced heights r,. The constant C determines the 
weighting between two pairs of centered differences. The 
limit of C = 1 corresponds to standard two-point finite differ- 
encing, which is generally accurate to second order. Further 
Taylor-series expansion, assuming a constant-spacing grid 
and a reasonably smooth function f(r), would give a result 
that is accurate to fourth order for C = 4/3. Note, though, that 
quantities in the ZEPHYR code are tabulated on a grid with 
variable spacing and they often exhibit numerical noise. In 
this case, errors may be reduced by essentially averaging be- 
tween the z'± 1 difference and the z'±2 difference (i.e., using 
C w 0.5). Tests with the ZEPHYR code found that C = 0.2 
provides the most accurate differentiation and noise reduc- 
tion, and this value is used in the models presented below. For 
the four zones at the bottom (; = 1 , 2) and top (; = N— 1 , N) of 
the grid we use a combination of the standard two-point finite 
difference expression and linear extrapolation to evaluate the 
derivatives. 

Figure 4 illustrates the convergence of the outer iteration 
process for the main polar coronal hole model described in 
§ 8.3. The displayed plasma parameters reach reasonably 
steady final values in about 70 outer iteration steps. The en- 
ergy convergence parameter (SE) decreases to its minimum 
range of variation (0.01-0.02) in about 100 iterations. Note, 
though, that (SE) is averaged over the entire radial grid; in 
many parts of the grid the convergence is much better than 
the average value indicates. Below the chromosphere-corona 
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transition region the convergence is excellent (i.e., (8E) is 
always less than 10~ 4 ). In the narrow transition region it- 
self, though, the piecewise smoothing discussed above smears 
out the temperature distribution so that the ideal "conduction 
dominated" solution cannot be sustained exactly. 9 Thus, the 
average value of the convergence parameter for 10 4 < T < 10 5 
K is about 0.2-0.3. In the corona and solar wind acceleration 
region above the transition region the average value of the 
convergence parameter is about 0.01, similar to the global av- 
erage. If the 30 or so grid zones of the transition region are 
excluded, the global average (SE) is reduced by about a factor 
of two to -0.005. 

The full ZEPHYR code comprises approximately 2000 
lines of Fortran. With the typical iteration parameters given 
above, the code runs in 2 to 4 hours of CPU time on various 
Sun Microsystems (Ultra and SunFire) computers. 

8. RESULTS 

In this section we present a series of solar atmosphere mod- 
els computed by the ZEPHYR code. A series of tests was first 
performed to evaluate the sensitivity of the dominant plasma 
properties to various input parameters (§§ 8.1-8.2). We used 
the grids of test models to determine the most likely input pa- 
rameters for a detailed model of a flux tube emerging from a 
polar coronal hole at solar minimum (§ 8.3). Then, using the 
same input parameters and varying only the radial dependence 
of the magnetic field Bo(r), we modeled the pole-to-equator 
variation of solar wind conditions at solar minimum (§ 8.4) 
and explored the properties of slow wind streams that are con- 
nected to flux tubes emerging from active regions (§ 8.5). 

In addition to the global magnetic field strength (our pri- 
mary "control knob") there are three key parameters that we 
varied in the course of exploring the physics of atmospheric 
heating and solar wind acceleration: 

1 . The photospheric acoustic flux Fsq injected at the lower 
boundary mainly affects the chromospheric heating. 
Probable values for Fsq seem to range between 10 7 and 
10 9 erg s -1 cm -2 (see, e.g., Musielak et al. 1994; Ulm- 
schneider et al. 1996, 2001; Carlsson & Stein 1997; 
Fawzy et al. 2002). For standard photospheric densi- 
ties and sound speeds, this gives a range for the photo- 
spheric acoustic velocity amplitude vmq of about 0.1 to 
1 km s" 1 . These values are a bit smaller than traditional 
"laminar" granulation velocities of 1 to 2 km s" 1 . The 
sources of propagating waves are believed to be con- 
centrated in the dark intergranular lanes and thus are 
able to extract only a fraction of the total kinetic energy 
of granulation (Rimmele et al. 1995; Nesis et al. 1997, 
1999; Cadavid et al. 2003). 

2. The photospheric Alfven wave amplitude v±q is spec- 
ified instead of the basal flux Fa or non-WKB ve- 
locity amplitude w±, since the latter quantities de- 
pend on the cancellation between upward and down- 
ward propagating waves that is determined as a part 
of the self-consistent solution. Observational determi- 
nations of w± from the footpoint motions of G-band 

9 Note that the solar atmosphere should also contain some degree of am- 
bipolar diffusion between ions and neutrals (Fontenla et al. 1990, 1991, 
1993). This effectively provides additional "smearing" of the otherwise ex- 
tremely sharp transition region. Our transition region thicknesses resemble 
those of Fontenla et al. by pure coincidence. None of the derived properties 
of the model either above or below the transition region appear to depend on 
this smearing. 
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FIG. 5. — Radial dependence of the background magnetic field Bo- The 
axisymmetric solar-minimum field of Banaszkiewicz et al. (1998) is shown 
for field lines originating at 9q = 0° (upper solid line), 16° (dashed line), 
24° (dot-dashed line), and 29.7° (dotted line). All have been modified at 
low heights using the model of Cranmer & van Ballegooijen (2005). Also 
shown are example active-region fields for /; = 0.03 and 0.07 i?© (dash-triple- 
dot lines) and the computed Alfven-wave magnetic amplitude B± for the 
polar coronal hole model (lower solid line). Above the largest height shown 
(r Ri 20 Rq) the magnetic field is nearly exactly radial, with Bo ^ r 2 - 

bright points yield mean values around 1 km s , with 
transient speeds up to 5 km s" 1 (e.g., Berger & Title 
1996). Larger-scale analyses of the horizontal diffusion 
of magnetic flux elements give smaller speeds of order 
0.1 to 0.3 km s" 1 (Schrijver et al. 1996). Recall that 
Vj_ is likely to be factors of 2 to 5 smaller than w±, so 
the range of possible values for Vx© may extend from 
below 0. 1 up to 1 or 2 km s" 1 . 

3. The photospheric Alfven wave correlation length L±@ 
sets the scale of the turbulent heating rate Q A (eq. |47l ). 
Once this parameter is set, the value of L± at all larger 
heights is determined by the adopted proportionality 

-1/2 

with B ' .A practical lower limit for L± Q seems to be 
about 10 km; i.e., the spatial scale over which radiative 
diffusion may inhibit the collapse of strong fields into 
thin flux tubes (e.g., Venkatakrishnan 1986; Sanchez 
Almeida 2001; Cameron & Galloway 2005). The up- 
per limit may be of the order of the size of photospheric 
granules (—1000 km). Possible intermediate length 
scales include the photospheric radius of a of thin flux 
tube (50-100 km), the width of an intergranular lane 
(about 300 km), and the mean separation between the 
flux tubes in the supergranular network (350-700 km; 
see Cranmer & van Ballegooijen 2005). 

All other parameters have been fixed with the values given in 
earlier sections. 

8.1. Coronal Parameter Study 

The test models discussed in this section all used the polar 
coronal hole magnetic field model that was derived by Cran- 
mer & van Ballegooijen (2005). In the photosphere, chromo- 
sphere, and low corona (i.e., from z = to 12 Mm), this model 
was obtained by tracing the radial magnetic field strength 
from the central axis of a two-dimensional numerical model 
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FIG. 6. — Contour plots of solar wind quantities resulting from varying 
the coronal heating parameters vj_q and Lj_q: ( a ) terminal wind speed Moo 
in units of km s _1 , (b) mass loss rate M in units of Mq yr -1 , (c) maximum 
coronal temperature in units of MK, (d) heliocentric critical radius in units of 
i?0. Also shown in each panel are the parameters chosen for the model of 
fast wind from a polar coronal hole discussed in § 8.3 (stars). SEE LAST 
PAGE OF PAPER FOR LARGER VERSION. 

of the supergranular network. This model contains thin inter- 
granular flux tubes between the photosphere and a "merging 
height" of 0.6 Mm where the flux tubes have expanded later- 
ally to the extent that the surrounding field-free plasma disap- 
pears. Above this height the merged network element under- 
goes further funnel-like horizontal expansion to fill a super- 
granular canopy. The field is directed completely vertically by 
the height of 12 Mm. Above this height we applied a slightly 
modified version of the empirically derived solar-minimum 
magnetic model of Banaszkiewicz et al. (1998). The radial 
dependence of Bq is shown in Figure 5. 

Figure 6 shows the result of producing a two-dimensional 
grid of ZEPHYR models by varying v^© and Lj_q. The 
acoustic flux F SQ was kept fixed at a median value of 10 8 
erg s" 1 cm" 2 (see above). The other two quantities were 
varied between the limits of 0.1 < v± Q < 1.5 km s" 1 and 
10 < L±q < 1000 km, with 9 points per quantity spread loga- 
rithmically between those limits. Figure 6a displays contours 
of the "terminal" outflow speed at the upper edge of the spa- 
tial grid, which we call u^. This quantity is slightly larger 
than the outflow speed at 1 AU, but never by more than 5%. 
Figure 6b shows the spherical mass loss rate, which is defined 
as 

M = Airpur 2 , (67) 

with the quantities on the right-hand side evaluated at the top 
of the grid (where A oc r 2 ). Figures 6c and 6d show the max- 
imum coronal temperature and the heliocentric radius of the 
wave-modified critical point (see § 6), respectively. 

Several general trends are evident in Figure 6. The mass 
loss rate is primarily dependent on vj_©, which represents the 
total amount of Alfvenic wave energy available to be dissi- 
pated in the corona. The outflow speed, though, seems to de- 
pend mainly on Lj_©; this parameter tells us where the wave 
energy is damped. For large values of L± & the damping oc- 
curs over a large range of heights, with increasingly more 
heating and acceleration taking place above the critical point. 
It has been known for some time that the relative heights of 
the critical point and the dominant energy deposition are key 
in determining the nature of a pressure-driven wind (Leer & 
Holzer 1980; Pneuman 1980; Leer et al. 1982). Heat that 
is deposited above the critical point is converted nearly com- 
pletely into kinetic energy of the wind (and a higher value 
of Mqo). On the other hand, low values of L± Q give a more 
concentrated heat deposition mostly below the critical point. 
This energy raises the temperature in the subsonic part of the 
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FIG. 7. — Coronal heating rates per unit mass (Qa/p) in units of erg s -1 g~' 
and temperatures (in K) for three models with constant values of L ± q = 75 
km and Fsq = 10 8 erg s _1 cirT 2 , and a range of values for vxq =0.1 (dotted 
lines), 0.255 (solid lines), and 0.65 km s -1 (dashed lines). 

corona, increases its scale height, and provides more down- 
ward heat conduction into the upper transition region. Less 
energy is thus available to accelerate the wind and 
lower. The mass loss rate M is generally believed to be set in 
the transition region by a balance between conduction, radia- 
tive losses, and an upward enthalpy flux (e.g., Hammer 1982; 
Leer et al. 1998). 

Interestingly, for a large range of parameters in the lower- 
right of Figure 6b (L± Q > 30 km and v± Q < 0.4 km/s), a 
combined power-law fit to the parameter dependence of the 
mass loss rate yields good agreement with the ZEPHYR mod- 
els for M oc v 2 [ 7 Q /L_ L 0. This resembles the classical Kol- 
mogorov (1941) energy flux due to isotropic hydrodynamic 
turbulence: v 3 /£. It seems to corroborate the general idea that 
the mass flux of a thermally driven wind is proportional to the 
deposited coronal energy flux. The large variation in coronal 
heating, though, does not seem to produce a large variation in 
the maximum coronal temperature. Figure 6c shows that this 
temperature varies only by about a factor of two over most 
of the parameter space. Owocki (2004) summarized how the 
energy losses due to both conduction and the solar wind act 
as an effective "thermostat" to keep the coronal temperature 
from varying too widely. 

Figure 7 illustrates how the coronal heating rate Q A changes 
as vj_ Q is varied and L± Q is kept fixed. We plot the heating 
rate per unit mass Qa/p in order to more easily show which 
heights receive the most heating on a particle-by-particle ba- 
sis. Note that larger values of vj_ produce larger values of 
both Q A I p and T and move their local maxima down to lower 
heights. The presence of damping, though, leads to differ- 
ences in how the heating rate depends on vj_ Q in various re- 
gions. In the lower corona (z « 0.05 Rq), before substantial 
Alfven-wave damping has had time to occur, the three models 
show a power-law dependence of Qa/p oc v^q that is similar 
to the mass loss rate dependence of M ex Vj/A for these mod- 
els. In the extended corona, though, the peak value of Qa/p 
(above z ~ 0.1 Rq) varies more weakly as v^q. The mod- 
els with higher Alfven wave amplitudes have undergone rela- 
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tively more damping in the extended corona than the models 
with lower v_l . 

The above arguments do not explain all of the features of 
the contours shown in Figure 6. For example, although there 
is a general trend (in the center and lower-right of Figure 6b) 
for a decreasing L±q to give a larger value of M, the glob- 
ally largest mass loss rate occurs in the upper-right corner of 
the plot — i.e., for the largest values of both v±q and L±q. 
In this region of the plot, added insight can be found in the 
"cold" wave-driven wind model of Holzer et al. (1983). For 
undamped Alfven waves that dominate the critical velocity 
(eq. |59l ) the terminal wind speed and mass loss rate can be 
computed analytically. Isolating the proportionality with v^ 
in this model yields 



M oc v 



X0 ' 



cx(l/v ; 



io)- 



constant . 



(68) 



where the constant term is typically negligible, thus giving 
oc v^'q (see eqs. [39] and [41] of Holzeret al. 1983). In the 
upper-right corners of Figures 6a and 6b, these relations come 
the closest to being satisfied. Specifically, for vj_q > 0.6 km 
s" 1 and L±q > 600 km, the best power-law fits to the model 
results are M oc v^g and v_l°q 9 . These parameters in- 

deed correspond to a wave-dominated critical velocity (i.e., 
large vj_ Q ) without much damping (large L±q). 

For the entire grid of values shown in Figure 6, a power- 
law is a poor fit to the overall variations of M and as 
a function of the two input parameters. For completeness, 
though, we report the best-fitting exponents over the full grid: 



M oc vY 7 L]°^ 7 and u c 



OC V 
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8.2. Chromo spheric Parameter Study 

We produced a second grid of exploratory ZEPHYR models 
that kept the coronal heating parameters (vj_0, Lj_©) fixed and 
varied only the basal acoustic wave flux Fsq. As described 
further in § 8.3, we chose optimal values of v^ Q = 0.255 km 
s" 1 and L±q = 75 km in order to model the fast wind that 
emerges from a polar coronal hole (see also the stars in Fig- 
ure 6). For the median value of Fsq = 10 8 erg s" 1 cm -2 used 
above, the transition region occurs at a relatively large height 
of z = 7200km w 0.01 Rq. A traditional view of chromo- 
spheric energy balance is that when the acoustic heating is 
increased, the height of the transition region should decrease 
(since the temperature would rise more rapidly as a function 
ofz). 

Figure 8 shows that the opposite actually occurs when we 
vary Fsq and keep everything else fixed. We produced a series 
of models with F SQ = 10 5 , 10 6 , 10 7 , 10 8 , 10 9 , and 10 10 erg s" 1 
cm" 2 , as well as a model with Fsq = (i.e., a model with 
Alfven wave heating only). The models with Fsq = and 10 5 
are virtually identical and the latter is not shown. 

Before discussing the chromospheric heating, we first note 
that the coronal and solar wind parameters for these models 
are remarkably constant. This shows that varying the acoustic 
heating has relatively little impact above the transition region. 
As Fsq is increased from zero to its maximum value, the mass 
loss rate decreases by only 7%, the terminal speed increases 
by less than 1%, and the peak coronal temperature decreases 
by 5%. 

Why does the transition region height increase as additional 
chromospheric heating is imposed? The models having larger 
values of Fsq have larger scale heights because they have 
both larger chromospheric temperatures and higher acoustic 
wave pressures. The larger scale heights lead to a more shal- 
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FIG. 8. — Temperature (left axis) and hydrogen number density (right 
axis) shown as a function of height above the photosphere (in km) for models 
that vary the acoustic wave flux Fsq (in units of erg s cm ) and keep the 
coronal heating parameters Vj and Lxq fixed (see captions). 

low density decrease in the approximately hydrostatic chro- 
mosphere. The transition region tends to occur at a critical 
density at which radiative cooling can no longer keep pace 
with the imposed acoustic heating (i.e., because A(T) has 
reached its peak and can increase no further). The "shallower" 
models with more acoustic heating reach this critical density 
at a larger height. 

The modeled transition region heights zjr range from 3800 
km (for F S q = 0) to 1 1500 km (for F SQ = 10 10 erg s" 1 cm" 2 ). 
For consistency, we define ztr as the height where T first 
reaches 2 x 10 5 K, which is the peak of the radiative cool- 
ing curve; see Figure 1. These heights are conspicuously 
larger than the standard values of, e.g., 1700 to 2300 km from 
"FAL" semi-empirical solar atmosphere models (Fontenla et 
al. 1993). Although there are some solar limb observations 
that suggest large values similar to our modeled zjr (Zhang et 
al. 1998), these may be affected by the dynamics of spicules 
and mass flows along closed magnetic loops, and thus may 
not be comparable directly to the present models (see also 
Filippov & Koutchmy 2000). 

The relative stretching of the chromospheres shown in Fig- 
ure 8 can be understood further by examining the density de- 
pendence of the total heating rate: (Qs+Qa) oc p v . The ex- 
ponent r\ was computed for all of the models by utilizing the 
radial dependence of both p and the total heating rate, with 



V 



dln(Q s + Q A )/dr 
dlnp/dr 



(69) 



Because the heating is balanced by radiative cooling (which 
usually has Qrad oc p 2 ), we find that values of r) closer to 2 give 
rise to a more extended "matching" between heating and cool- 
ing at chromospheric temperatures — and thus a larger overall 
extent of the chromosphere. Thus, since rj < 2 for all of the 
models shown below, we can understand how larger [smaller] 
values of T) correspond to a higher [lower] transition region 
height. For the shock heating described in § 4.1, the possible 
values of 77 range between -0.5 and +1 . The lower limit would 
occur for an undamped weak shock train, for which T AS oc vl 
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(eq. ll33lD and pv^ = constant. The upper limit occurs when 
the shock train strengthens to the point where vm saturates to 
a nearly constant value and thus TAS « constant. (The peak 
Mach number M\ tends to saturate at values between 2 and 3 
for the majority of the Fsq models.) 

The ZEPHYR models shown in Figure 8 exhibit interme- 
diate values of the density exponent 77, with the value of the 
exponent varying slightly as a function of height in the chro- 
mosphere. The lowest value of z TR corresponds to the F S q = 
model and 77 = 0.4-0.5. The highest value of z TR corresponds 
to the Fsq = 10 1() model and 77 = 0.8-1 . It is clear that, in these 
models, an even lower (FAL-like) transition region height 
would require r\ < 0, and thus would require the heating rate 
to increase with height. Such heating may occur for the fol- 
lowing reasons. 

1 . As described above, if the acoustic wave energy were 
to act more like a weak undamped shock train, r\ would 
approach -0.5. The emerging picture of the chromo- 
sphere as a confluence of multiple acoustic sources 
(each expanding quasi-spherically from small regions 
at or below the photosphere) may be consistent with 
the need for more low-amplitude acoustic wave energy 
(see, e.g., Ulmschneider et al. 2005). A collection of 
weak and incoherent acoustic wave packets may also 
exhibit effectively larger frequencies because of their 
random phases and constructive interference. 

2. There could be additional important sources of linear 
wave dissipation in the low chromosphere, similar in 
form to 7cond (which gives i)«-l to -0.5 depending on 
the magnetic geometry). For example, Goodman (2000, 
2004) suggested that the resistive dissipation of currents 
driven by MHD waves may provide enough energy to 
heat the chromosphere. 

3. Similarly, it has been argued on the basis of multiple 
lines of observational evidence that both solar and stel- 
lar chromospheres may be dominated by magnetic and 
not acoustic heating (e.g., Judge & Carpenter 1998; 
Judge et al. 2003; Bercik et al. 2005). If these mech- 
anisms grow stronger as one moves from the lower to 
the upper chromosphere (the latter being more strongly 
ionized and magnetized), this could produce the radial 
increase of the total heating rate that would be needed 
to push down the transition region. 

However, we should note that Anderson & Athay (1989a) 
determined empirically that the chromospheric heating rate 
seems to be roughly proportional to p (and thus r\ w 1) 
throughout most of the chromosphere. The mechanisms in- 
cluded in the ZEPHYR models may thus be reasonable with- 
out the need for additional physics. 

We experimented with several other variations on the 
acoustic heating in order to test the assumptions described 
in § 4. A test model was created with the power spectrum 
Ps(u>) divided into a finer mesh of discrete frequency bins: 
separated by factors of V2 rather than 2. This model, which 
had Fsq = 10 8 erg s" 1 cm" 2 , produced a nearly identical chro- 
mosphere to the standard model with the same value of Fsq. 
The fine-mesh model's transition region height zjr was about 
10% higher than that of the standard model because slightly 
more power was given to the lower frequencies (which form 
shocks at larger heights). Presumably, using a coarser fre- 
quency spectrum would lead to a lower value of ztr, but 



this would start to become unfaithful to the modeled spec- 
tral shape (eq. l42l ). Another model was run with just a sin- 
gle frequency bin; i.e., a monochromatic low-frequency wave 
train with u> = ui- dc . This model was designed to explore what 
would happen if the high-frequency tail in Ps(ll>) were absent 
(e.g., Fossum & Carlsson 2005, 2006). This model exhibited 
Zjr = 14300 km, about a factor of two larger than the standard 
model, and would certainly require some additional kind of 
chromospheric heating to produce a realistically low value of 
Ztr- 

Note that for all of the above models the minimum tem- 
perature in the upper photosphere never dropped below the 
radiative equilibrium value of r ra d w 4500 K. There is obser- 
vational evidence, though, for a lower minimum temperature 
that could extend intermittently down to values between 3000 
to 4000 K (see model A of Fontenla et al. 1993; Carlsson 
& Stein 1997; as well as recent work by Ayres et al. 2006; 
Fontenla et al. 2006). Additional sources of cooling that could 
be included in ZEPHYR include molecular opacity and dust 
formation (for cooler stars) as well as adiabatic expansion ef- 
fects due to waves and shocks that may not be confined to the 
modeled flux tube. 

8.3. Polar Coronal Hole Model 

The exploratory models discussed above led to a choice 
for the optimal set of parameters which would reproduce the 
observed properties of high-speed solar wind streams that 
emerge from polar coronal holes (mainly at the minimum of 
the solar cycle). These parameters are v^© = 0.255 km s -1 , 
L±q = 75 km, and Fsq = 10 8 erg s -1 cm" 2 . The terminal speed 
Woo and mass loss rate M computed for this model are 753.5 
km s" 1 and 1.88 x 10~ 14 MQyr _1 , respectively. These values 
give a total hydrogen number density of 2.9 cm" 3 at r = 1 AU. 

The semi-empirical Alfvenic turbulence model of Cranmer 
& van Ballegooijen (2005) determined a value for L± at the 
mid-chromosphere merging height of about 1100 km, which 
corresponds to a photospheric value L±q w 320 km. The 
Cranmer & van Ballegooijen (2005) model also predicted a 
photospheric non-WKB amplitude w±q of about 3.1 km s -1 , 
which translates into vj_ 0.46 km s . These values are 
both slightly higher than the ones chosen on the basis of the 
self-consistent ZEPHYR models, though they certainly fall 
within the range of plausibility. The smaller value of the tur- 
bulence correlation length (75 km) seems to be consistent with 
the observed horizontal size of a thin flux tube in the photo- 
sphere (50-100 km). This seems nicely consistent with being 
an outer "stirring" scale for the turbulence, since the horizon- 
tal shaking and distortion of the flux tube can be expected to 
take place mainly on the spatial scale of its own size. (Ear- 
lier justifications for the larger [~ 300 km] correlation length 
were based on this being either the mean separation between 
flux tubes or the intergranular lane width. These scales may 
be excited by convective driving as well, but it makes sense 
for the primary response of each tube to be on the smaller 
scale of its radius or diameter.) 

For simplicity we continue to use the fiducial value for Fsq 
that was used in the two-dimensional grid shown in Figure 6. 
This energy flux density (10 s erg s" 1 cm" 2 ) corresponds to a 
photospheric acoustic wave amplitude v\\q of about 0.29 km 
s -1 , and it is close to that computed by Musielak et al. (2000) 
for linear longitudinal flux-tube waves. The flux computed 
in the Lighthill-Stein sound generation models depends sen- 
sitively on the partitioning of gas and magnetic pressure in 
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FIG. 9. — Radial dependence of various plasma parameters for a polar coro- 
nal hole, (a) Velocity quantities in km s : wind speed u (solid line), Alfven 
speed V A (dashed line), Alfven wave amplitude Vx (dotted line), acoustic 
wave amplitude (vii) (dot-dashed line), and the wave-modified critical point 
(diamond symbol), (b) Temperature T in MK (solid line), spectrum-averaged 
reflection coefficient TZ (dashed line), and mass density p (dot-dashed line; 
right axis), (c) Areas denote the relative contribution of terms to the energy 
conservation equation; see labels. (Q at i v denotes the advection terms on the 
left-hand side of equation ^3), other terms are defined in the text.) 

the photosphere; these models have Fsq oc (Bo/B eq )~ p , where 
B eq = (87rP) I//2 is the field strength consistent with equipar- 
tition between gas and magnetic pressures, and p takes on 
large values typically between 6 and 9. Our polar coronal 
hole model has a photospheric ratio Bq/B^ = 0.77. Interpo- 
lating from Table 1 of Musielak et al. (2000) would yield a 
flux for this model of about 1.1 x 10 8 erg s -1 cm" 2 , which is 
extremely close to what we use. 

Figure 9 displays the radial dependence of several key 
plasma parameters for the polar coronal hole model. The 
critical point is denoted by a symbol, with r c = 1.84^ Q and 
u c = 166 km s" 1 . The wave-modified critical point is only 
a small distance above the classical sonic point, which oc- 
curs at 1.74 Rq (where u = a = 149 km s" 1 ). This is some- 
what surprising because the fast solar wind is typically as- 
sumed to be strongly "wave dominated," and thus one might 
have expected the modified critical point to be far from the 
unmodified sonic point. Both of these points are also near 
the so-called turbopause radius at which the bulk solar wind 



speed begins to exceed the turbulent fluctuation amplitude v± 
(r = \ .5ARq\ see Veselovsky 2001). In contrast, the Alfvenic 
singular point (where u = Va) occurs at a substantially larger 
radius of ta = 10. 8Rq and ua = 509 km s" 1 . Note that Figure 
9a displays a frequency-integrated acoustic wave amplitude 
that we define for convenience as 



1/2 



V 



(70) 



bins 



(see eq. 12411 ). where we summed over the discrete frequency 
bins discussed in § 4.2 and assumed s = 2 everywhere. This 
acoustic wave amplitude peaks at about 17 km s" 1 in the upper 
chromosphere (and is mildly nonlinear, with (y\\)/c s peaking 
there at 0.96), then is damped strongly in the low corona. 

The temperature and density curves shown in Figure 9b ap- 
pear to be in reasonable agreement with prior expectations for 
the fast solar wind (e.g., Kohl et al. 2006). When p is con- 
verted into electron number density n e , the modeled values 
agree with the canonical polar coronal hole measurements of 
Sittler & Guhathakurta (1999) to within ± 20%. These mea- 
surements come mainly from white-light polarization bright- 
ness (pB) observations between r = 1.5 and 4 Rq. Consider- 
ing that the absolute pB calibration uncertainties are typically 
about 10% and that the measured density depends on the fill- 
ing factor of polar plumes along the line of sight (which in- 
troduces pB variations up to a factor of two in magnitude; see 
Fisher & Guhathakurta 1995), this agreement is good. 

The Alfven wave amplitudes (vj_ in Figure 9a and B± in 
Figure 5) compare favorably with the semi-empirical mod- 
els of Cranmer & van Ballegooijen (2005), who compared 
those models with various remote-sensing and in situ mea- 
surements. The spectrum-averaged reflection coefficient 1Z 
shown in Figure 9b is also similar to that of Cranmer & van 
Ballegooijen (2005); note that it need not be monotonically 
decreasing with increasing distance. The photospheric value 
of TZ in this model is 0.883, which is smaller than the value 
of 0.974 from Cranmer & van Ballegooijen (2005). The latter 
model utilized the FAL-C' temperature structure which had a 
sharper transition region than that produced by the ZEPHYR 
code. A more abrupt transition region has a larger Va gra- 
dient and thus experiences stronger reflection. Defining the 
transition region thickness <5ztr as the distance between tem- 
peratures of 2 x 10 4 and 2 x 10 5 K, the FAL-C model has 
fern =120 km and the ZEPHYR polar coronal hole model 
has a value of 560 km. Although some of the computed thick- 
ness may be due to numerical smoothing (see § 7) not all of it 
is. The ZEPHYR code did produce a sharper transition region 
(fern = 270 km) for the Fsq = model shown in Figure 8 — 
which was also the model with the lowest transition region 
height z TR . 

The radial dependence of temperature in the extended 
corona and heliosphere can be compared with various ana- 
lytic limiting cases (e.g., Hundhausen 1972). At large dis- 
tances where the wind speed is approximately constant and 
p oc r" 2 , there are several cases that give a power-law depen- 
dence T(r) oc r~* '. The polar coronal hole model shown in 
Figure 9 exhibits (3 ~ 0.30 at 1 AU. If the advection terms 
on the left-hand side of equation (O are dominant, one ob- 
tains the purely adiabatic exponent (3 = 4/3. If classical heat 
conduction were to dominate the energy equation out to 1 
AU, (3 = 2/1. However, we find that conduction and advec- 
tion tend to balance one another in the heliosphere and yield 
intermediate values. For the case that the advection terms 
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are balanced by classical Spitzer-Harm conduction, there is 
no exact power-law solution for T(r). For advection being 
balanced by collisionless conduction (eq. l22l ). one obtains 
(3 = 8/(6 + 3a c ), which for the adopted value a c = 4 gives 
(3 = 4/9. In the ZEPHYR models, the hybrid heat conduction 
(eq. fl"4l ) is mostly collisionless at 1 AU, but it still contains 
some classical heat conduction. The computed values of j3 
thus tend to fall between the two above cases of 2/7 and 4/9. 

Figure 9c shows the terms that dominate the energy con- 
servation equation as a function of distance. The areas plotted 
here were computed by normalizing the absolute values of the 
individual terms in equation (01 by the maximum value at each 
height (as in the denominator of eq. 11631 ) and then "stacking" 
them so that together they fill the region between and 1. 
The photosphere is essentially definable as the region where 
radiative heating balances radiative cooling (both denoted as 
2rad) and the chromosphere is the region where acoustic heat- 
ing (Qs) balances £> ra d- The transition region and the very low 
corona (T < 0.5 MK) exhibit a complicated balance of ra- 
diation, conduction, advection (i.e., enthalpy flux), and some 
Alfvenic and acoustic heating. The extended corona is mainly 
a balance between the Alfven wave heating {Qa) and the ad- 
vection terms on the left-hand side of equation ||3}, although 
conduction remains nonnegligible. In the heliosphere above 
20-30 Rq the direct heating becomes less important and ad- 
vection balances conduction. 

In order to produce further comparisons with in situ solar 
wind measurements, we computed the nonequilibrium ion- 
ization balance of oxygen as a function of distance in the 
ZEPHYR models. Specifically, the ratio of number densities 
of 7+ to O s+ is often used to aid in the identification of fast 
and slow wind streams (e.g., Zurbuchen et al. 2002). Because 
of the steep decline in electron density with increasing height, 
solar wind ions above a certain "freezing-in radius" encounter 
virtually no electrons, and thus are not sensitive to ionization 
and recombination processes in interplanetary space (Hund- 
hausenetal. 1968; Owocki et al. 1983). Interplanetary charge 
states thus carry information about the plasma properties in 
the corona. We adopted the nonequilibrium ionization code 
of Gaetz et al. (1988), Esser et al. (1998), and Esser & Edgar 
(2000, 2001) to the ZEPHYR models at temperatures above 
10 4 K. Ionization and recombination are most sensitive to the 
electron velocity distribution in the corona, and we relied on 
our basic one-fluid assumption of a Maxwellian distribution 
with T = T e = T p . We also assumed that all oxygen ions flow 
with the bulk wind speed m, independent of their charge. 

Figure 10 shows the result of computing the oxygen ioniza- 
tion state for the polar coronal hole model. The mean charge 
state (Zo) is computed as an average of the net charge of the 
ions (in units of e) weighted by the computed number density 
fractions of each stage of ionization. The equilibrium solution 
assumes a local "coronal" balance between collisional ioniza- 
tion, radiative recombination, and dielectronic recombination, 
and is a strict function of temperature. The full nonequilib- 
rium solution freezes in at a relatively low height in the corona 
(z ~ 0.05 Rq) and is roughly constant at all larger heights. 
Note, though, that the equilibrium and nonequilibrium solu- 
tions are also different from one another in the upper transi- 
tion region (z « 0.01-0.02 Rq) where the wind flow time over 
a scale height is beginning to approach the relevant ionization 
and recombination times. The dominant ionization state at 1 
AU is 6+ (i.e., (Z ) w 6) and the ratio of 7+ to O e+ is dis- 
cussed further below. 
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FIG. 10. — Oxygen ionization versus radial distance for the polar coronal 
hole model. The mean charge states for the nonequilibrium "frozen in" model 
(solid line) and the local coronal equilibrium (dashed line) are compared with 
the logarithm of the temperature T, in K (dotted line). 

8.4. Solar Minimum Axisymmetric Field 

There is a definite empirical relationship between the solar 
wind speed measured in situ and the inferred lateral expan- 
sion of magnetic flux tubes near the Sun (Levine et al. 1977; 
Wang & Sheeley 1990; Arge & Pizzo 2000). Specifically, flux 
tubes that expand more rapidly between the coronal base and 
r 2.5 Rq tend to have lower wind speeds at 1 AU. There 
have been many theoretical attempts to explain this observed 
antic orrelation. Some models have treated the flux expansion 
as the fundamental cause of the difference between fast and 
slow streams (e.g., Kovalenko 1978, 1981; Wang & Sheeley 
1991, 2006; Wang 1993; Cranmer 2005a), some have treated 
it as a by-product of varying boundary conditions in the low 
corona (e.g., Fisk 2003; Schwadron & McComas 2003), and 
some have considered a combination of the two (Bravo & 
Stewart 1997; Suzuki 2006). 

Our goal is to test these ideas by varying only the flux 
expansion rate and keeping the lower boundary parameters 
(vj_0, L±q, Fsq) fixed. The first way that we modify the flux 
expansion is to utilize the full two-dimensional axisymmet- 
ric magnetic field model of Banaszkiewicz et al. (1998). This 
model utilizes a sum of dipole, quadrupole, and current-sheet 
(effective monopole) terms to reproduce various observed 
properties of the solar-minimum field. Note that other as- 
sumptions lead to slightly different solar-minimum field con- 
figurations (e.g., Sittler & Guhathakurta 1999, 2002; Vasquez 
et al. 2003) but the general trends of the polar and equato- 
rial expansion factors are not substantially different from the 
Banaszkiewicz et al. (1998) case. 

Figure 1 1 illustrates a selection of the open field lines in the 
Banaszkiewicz et al. (1998) model. Figure 5 also shows the 
magnetic field Bo(r) traced along a subset of these field lines. 
We define the surface colatitude 9q as the main identifier of 
each flux tube. The polar flux tube that was used above in 
§ 8.3 has do = 0. We take the "last" open field line (at the 
outermost edge of the open-field region) as 9q = 29.7° '. Be- 
cause the latter field line eventually stretches to a colatitude 
of nearly 90° we also call this the "equatorial" flux tube. We 
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FIG. 11. — Upper panel: axisymmetric magnetic field geometry of 
Banaszkiewicz et al. (1998) with the radii of wave-modified critical points 
marked by open diamond symbols. Lower panel: outflow speeds for a series 
of axisymmetric flux tube models (surface colatitudes 6q given in captions) 
with critical point heights denoted by open diamonds, and the Alfven wave 
amplitude v± for the polar (p) and equatorial (e) models (dotted lines). 

assume that the differences in Bo between flux tubes occur 
only above z ~ 0.01 Rq and that in the photosphere and chro- 
mosphere all flux tubes are identical (see, e.g., Aiouaz & Rast 
2006). 

Let us define an effective Wang & Sheeley (1990) superra- 
dial expansion factor for each flux tube as 



/ss 



(fi()'" 2 )base 

CBor 2 ) ss 



(71) 



where the coronal base is fixed at rb ase = 1 -04 Rq (i.e., a height 
equivalent to the size of one supergranule, well above any 
"canopy" structure) and we use the traditional source surface 
radius of r ss = 2.5Rq (Hoeksema & Scherrer 1986). It would 
not be proper to use the photosphere as the basal height be- 
cause the low-resolution magnetograms that typically are used 
to determine / ss resolve neither the individual photospheric 
flux tubes nor the contrast between supergranular network and 
cell interior. For the Banaszkiewicz et al. (1998) field configu- 
ration, / ss ranges from 4.5 for the polar flux tube to 9. 1 for the 
equatorial flux tube. Note that this model does not produce 
the larger values of / ss 20-40 that have been found during 
more active phases of the solar cycle (see § 8.5). 

A series of 20 ZEPHYR models was computed with the 
same photospheric boundary conditions as the polar coronal 
hole model discussed in § 8.3, but varying do from to 29.7°. 
Going from pole to equator, the asymptotic solar wind speed 
Uoo was found to decrease and the mass loss rate M was found 
to increase. Figure 1 1 shows how the radial dependence of 
outflow speed changes as Oo increases from to 29.7°. There 
is an abrupt bifurcation between two classes of solar wind 
solution: (1) polar solutions that have fast wind speeds, low 
densities, and wave-modified critical points r c < 2R Q , and (2) 
equatorial solutions that have slow wind speeds, high densi- 
ties, and r c > 4-Rq. The colatitude at which this bifurcation 



occurs is 0q ~ 24.25°. As discussed by Vasquez et al. (2003) 
and in § 6 above, when there is more than one possible loca- 
tion for the critical point, the most stable time-steady solution 
tends to be the one with the largest value of r c . For the Ba- 
naszkiewicz et al. (1998) field geometry, the second (large-r c ) 
solution appears only for Oq > 24.25°. The field line that di- 
vides the two classes of solutions extends out to a heliospheric 
latitude (measured from the ecliptic plane at r ^> Rq) of 16°. 

Figure 1 1 shows that the solar wind solutions all exhibit 
substantial wind speeds just above the transition region (u « 
10 km s" 1 ). Esser et al. (2005) found that such rapid flows 
at the coronal base are a consequence of the rapid superra- 
dial divergence in supergranular funnels, and that these flows 
seem to be necessary to explain observed H I Lya emission 
from the solar disk. Note also that our modeled flux tubes 
nearest to the equator exhibit a pronounced local minimum in 
the wind speed at the cusp-radius of the current sheet in the 
Banaszkiewicz et al. (1998) model. Observational evidence 
for this kind of "stagnation point" in equatorial streamers was 
reported by Strachan et al. (2002) and was also modeled in a 
similar way as in this paper by Chen et al. (2004) and Li et al. 
(2004). 

It is also evident from Figure 1 1 that the polar model has 
a larger Alfven wave amplitude than the equatorial model. 
Some of this difference can be attributed directly to the vary- 
ing flux-tube area factors A(r) (i.e., different manifestations of 
wave action conservation). However, the equatorial flux tube 
also exhibits relatively more damping than the polar model. In 
the corona, the Alfven waves essentially "spend more time" 
at low heights in the equatorial flux tube because of its lower 
phase speed (u + Va), and this allows a given damping rate 
Qa to have more of an impact on diminishing the wave en- 
ergy. Thus, the equatorial models exhibit more heating at low 
heights (due to the increased damping and the higher densi- 
ties at the coronal base) and less heating at large heights (due 
to the lower wave amplitudes vj_ and larger turbulent correla- 
tion lengths L± in the extended corona). The peak tempera- 
ture for the polar model, T = 1 .35 MK, occurs at r = 1 .86/?0, 
and the peak temperature for the equatorial model, T = 1.29 
MK, has a similar magnitude but occurs at a lower height of 
r= 1.267?q (see also Figure 17below). This shift in the range 
of heights over which coronal heating occurs seems to be an 
important factor in producing a fast solar wind for the polar 
model (which is heated in the supersonic region) and a slow 
wind for the equatorial model (which is heated in the subsonic 
region); see § 8.1. The differences in temperature persist out 
to 1 AU, where the polar temperature (0.41 MK) exceeds the 
equatorial temperature (0.16 MK) by a larger relative amount 
than in the extended corona. The radial temperature exponent 
j3 (i.e., T cx r' 13 ) grows steeper from pole (f3 « 0.30) to equa- 
tor (/? w 0.47). 

Although the photospheric Alfven wave amplitude param- 
eter v^o was kept fixed in these models, the degree of non- 
WKB reflection decreased slightly from the polar to the equa- 
torial flux tube. This resulted in a decrease in the computed 
value of w±q (eq. Il56l ) from its polar value of 0.72 km s" 1 to 
an equatorial value of 0.59 km s" 1 . Interestingly, if we could 
have forced w±q to remain fixed as a function of latitude, the 
resulting equatorial value of vj_ would have been about 20% 
larger than the polar value. Figure 6 indicates that the effect 
would have been an even larger equatorial mass loss rate and a 
slightly lower u^. Thus, the current set of models is in some 
sense "robust" because these effects arise even without such 
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FIG. 12. — Latitudinal dependence of (a) outflow speed, (b) number den- 
sity, and (c) temperature measured in interplanetary space (see text for de- 
tails). Data from the Ulysses first polar pass (thin solid lines) are compared 
to the standard ZEPHYR models computed along axisymmetric superradial 
flux tubes (thick solid lines) and the 'Durham' ZEPHYR models computed 
along the same flux tubes (thick dashed lines). 

an extra fine-tuning in w±q. 

The detailed latitudinal dependence of the ZEPHYR mod- 
els (for an observer in interplanetary space) is shown in Fig- 
ure 12. We compare the model predictions with Ulysses 
SWOOPS (Solar Wind Observations Over the Poles of the 
Sun; Bame et al. 1992) measurements during its first fast lat- 
itude scan in 1994 and 1995 (see Goldstein et al. 1996). The 
measured values for flow speed, proton number density, and 
proton temperature were obtained at heliocentric distances be- 
tween 1.3 and 2.3 AU. The densities were scaled to a com- 
mon distance of 1 AU by using an assumed r~ 2 dependence. 
For comparison we plot Uoo, which does not vary substan- 
tially between these distances, and the proton number density 
(n p = p/oth/1.2) at 1 AU, scaled for an assumed helium-to- 
hydrogen number density ratio of 0.05. The model tempera- 
ture T(r) has been interpolated to varying heights as a function 
of latitude to match the elliptical orbit of the Ulysses probe 
during its polar pass. Also, we used the Banaszkiewicz et al. 



(1998) model to map the surface colatitude 9q of each model 
flux tube to its corresponding heliospheric latitude at r — > oo. 
The latter is the abscissa coordinate used in Figure 12. To 
within about 10% accuracy, one can estimate this latitude as 
|9O°-30 O |. 

There are two sets of pole-to-equator ZEPHYR models 
shown in Figure 12. We compare the standard model dis- 
cussed above with a preliminary model presented by Cran- 
mer (2006), which we hereafter call the "Durham" model. 10 
The models differ in two specific ways. (1) The standard 
model used the Alfven wave frequency spectrum shown in 
Figure 3, but the Durham model assumed a single frequency 
for the Alfven waves — corresponding to a period of 6 min- 
utes. (2) The standard model used an exponent n = 1 in 
equation (l48l . but the Durham model used n = 2. The lat- 
ter was an early guess for the functional form of £ tul -b that 
we used prior to the incorporation of the insights of Oughton 
et al. (2006). The larger exponent resulted in more efficient 
quenching of the MHD turbulence in the lower atmosphere, 
which thus required a larger amount of Alfven wave energy 
to produce plasma conditions appropriate for a polar coro- 
nal hole. Specifically, the Durham model used vj_ Q = 0.42 
km s" 1 , L±q = 120 km, and the same value of F$q as was 
used in the standard model. The Durham model's and M 
for the polar flux tube were nearly identical to those of the 
standard model. The peak coronal temperature, though, was 
substantially larger (2.06 MK) and the temperature at 1 AU 
was slightly smaller (0.31 MK) in comparison to the standard 
polar model. 

There is reasonably good agreement between the modeled 
and observed plasma properties shown in Figure 12. We 
of course do not reproduce any of the measured north-south 
asymmetry in the data, except for that in T p that arose because 
the path of Ulysses was not perfectly north-south symmetric. 
We also do not attempt to model the rapid near-ecliptic vari- 
ability that seems to result from longitudinal variations in Bq. 
In some sense, the Durham model does better at achieving the 
low equatorial wind speeds (down to 330 km s" 1 ) and high 
number densities (up to about 10 cm" 3 ) that are often ob- 
served in the plane of the ecliptic, but that model exhibits a 
narrower range of slow-wind (large-r c ) latitudes than both the 
standard model and the Ulysses data. It is expected, though, 
that some fraction of the observed ±20° latitudinal extent of 
the slow-wind region may be due to a slight tilt in the main 
dipole component of the solar magnetic field (i.e., a "balle- 
rina skirt" phenomenon) because the first Ulysses polar pass 
did not occur exactly at the minimum of the activity cycle. 

Note from Figure 12c that the ZEPHYR model tempera- 
tures between 1.3 and 2.3 AU are systematically higher than 
the measured values of T p . The trend as a function of the 
coupled distance-latitude coordinate, though, is similar to the 
measurements. Our heliospheric temperatures are sensitive to 
the adopted method of bridging between collisional and col- 
lisionless heat conduction (eq. lfT4l ) and the present method 
is admittedly somewhat ad hoc. However, some of the dis- 
crepancy in Figure 12c may come from the fact that T p ^T e 
in interplanetary space and that we are essentially modeling 
their average. In situ measurements have shown that T p > T e 
in the fast solar wind and T p < T e in the slow wind (see, e.g., 
Feldman & Marsch 1997). The discrepancy is worst in the 

10 This model was presented at the 37th meeting of the Solar Physics Di- 
vision of the American Astronomical Society in Durham, New Hampshire 
(25-30 June 2006). 



24 



CRANMER, VAN BALLEGOOIJEN, & EDGAR 



800 



700 



600 



500 



400 



300 





i i i i i i 

Bana. (standard) - 




1 Bana. (Durham) : 

i active B ■ 




1 \ Y>. : 








Wang & Sheeley (1990) " 
, , i , , 



7.5 - 



10 



50 



FIG. 13. — Terminal outflow speeds for the various se ts o f ZEPHYR mod- 
els plotted versus the flux-tube expansion factor / ss (eq. 1711 ). compared with 
the empirically derived relationship from Table 2 of Wang & Sheeley (1990). 
The standard (solid line) and Durham (dashed line) pole-to-equator models 
computed using the Banaszkiewicz magnetic field are shown along with the 
active-region extensions to equatorial Bq(t) discussed in § 8.5 (dot-dashed 
lines). The polar (open circle) and equatorial (filled circles) axisymmetric 
models are highlighted. 

polar, fast-wind regions where our modeled average temper- 
ature would exceed both the observed T p and T e . However, 
the ZEPHYR code does not yet contain the dissipative physi- 
cal processes that would lead to temperature equilibration be- 
tween protons and electrons; these may alter subtly both the 
net heating and the form of the conduction flux qu. 

Figure 13 plots the terminal speeds of the standard and 
Durham ZEPHYR models as a function of the Wang-Sheeley 
superradial expansion factor / ss . The abrupt changes in Uoo 
that occur because of the emergence of the outer critical point 
are evident as vertical discontinuities at the critical values of 
9(,. Because of the relatively narrow range of / ss values sub- 
tended by the pole-to-equator (Banaszkiewicz et al. 1998) flux 
tubes, we see that the standard model actually follows the the 
empirical Wang & Sheeley (1990) relationship better than the 
Durham model, despite the fact that the standard model has a 
relatively high equatorial terminal speed of 500 km s -1 . 

Figure 14 compares modeled and measured Alfven wave 
amplitudes in the heliosphere as a function of the in situ wind 
speed. It is difficult to process the richly nonlinear turbulent 
fluctuations (as measured in situ) into a single "amplitude," 
so we made an attempt to convert the ZEPHYR values of 
vj_ into something analogous to the frequency-averaged El- 
sasser spectrum quantities reported by Tu et al. (1992). Those 
data were obtained by the two Helios spacecraft between 0.29 
and 0.52 AU in the years 1979-1980. (Measurements taken 
closest to the Sun are optimal if the goal is to obtain the 
properties of the "pristine" Alfven waves.) Individual time 
series of Z^/2 (outward-propagating kinetic energy) magni- 
tudes for one-day periods were transformed into power spec- 
tra, denoted e + (f). The values shown in Figure 14 were av- 
eraged over frequencies / between 10~ 4 and 2 x 10~ 4 Hz. 
The plotted data points also correspond to time periods when 
the solar wind exhibited its highest range of total energy flux 
(Plot = pu[u 2 + V esc]/2X which we verified to correspond most 
closely with the ZEPHYR model results. 

We extracted the Elsasser amplitudes Z_ from the models 
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FIG. 14. — Comparison of Helios Elsasser spectra e + at r = 0.3-0.5 AU 
with ZEPHYR model predictions at a mean distance of r = 0.4 AU. Data 
from Figure 3 of Tu et al. (1992) is shown (gray points). Standard (solid 
line) and Durham (dashed line) pole-to-equator models are plotted together 
with active-region flux tube models (dot-dashed lines) in a similar way as in 
Figure 13. 

at a mean distance of r = 0.4 AU and converted them into the 
average value of <? + as follows. First, representative spectral 
shapes for these fluctuations were used in order to determine 
how much of the total energy is contributed by the specified 
frequency band (1-2 x 10~ 4 Hz). We used the power-law ex- 
ponents and breakpoints reported by Tu et al. (1989), using 
the same Helios data, to determine that fast and slow wind 
streams contain approximately 11% and 21%, respectively, 
of their power in this band. We used a mean value of 16% 
for all ZEPHYR models to avoid introducing any possibly 
spurious correlations with wind speed. To convert the (im- 
plicitly frequency-integrated) energy fraction into an "aver- 
age" spectrum quantity we divided by an effective frequency 
/ e ff that was defined to perform this conversion exactly for 
the idealized spectra of Tu et al. (1989); it turned out that 
/eff ~ 9.5 x 10~ 5 Hz for both fast and slow wind spectra and 
we used just 10~ 4 Hz. Thus we plot the derived quantity 

+ 0.16(Z?/2) 



KHHz 



kmV 2 Hz _ 



(72) 



in Figure 14 in order to compare with the measured average 
spectra from Tu et al. (1992). The modeled and measured 
values — and trends with wind speed — appear to agree quite 
nicely. 

In the remainder of this section we discuss two supplemen- 
tary calculations of in situ quantities that are often used to 
diagnose the properties of fast and slow wind streams: the 
ionization fraction (Figure 15) and the FIP effect (Figure 16). 

1 . The nonequilibrium oxygen ionization state at 1 AU has 
been computed for each of the standard and Durham model 
flux tubes. Figure 15 shows the ratio of 7+ to O e+ number 
densities as a function of the modeled terminal speed. These 
are compared to statistical summaries of Ulysses SWICS (So- 
lar Wind Ion Composition Spectrometer; Gloeckler et al. 
1992) measurements of this ratio taken during two time pe- 
riods: the first polar pass discussed above (September 1994 
to August 1995) and the initial phase of the mission that sam- 
pled the peak of the solar cycle (December 1990 to September 
1994). The data points, obtained from the online database of 
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FIG. 15. — Oxygen freezing-in ratio (i.e., ratio of O to O number 
densities) at 1 AU plotted as a function of solar wind speed in interplanetary 
space. Symbols and line styles for the models are the same as in Figure 
14. Binned Ulysses data from the 1990-1994 solar maximum time period 
(light gray region) and the 1994-1995 fast latitude scan (dark gray region) 
are shown for comparison. 

T. Zurbuchen and R. von Steiger (see also von Steiger et al. 

2000) , were grouped in sequential bins of wind speed of size 
25 km s" 1 , and the mean and standard deviation of the ioniza- 
tion ratio in each bin were computed. The regions plotted in 
Figure 15 show the values within ±1 standard deviations of 
the mean for each bin in the two data sets. 

The general trend of the pole-to-equator ZEPHYR models 
is for lower wind speeds to have higher ratios of 7+ to 6+ 
at 1 AU, in qualitative agreement with the observations. The 
absolute values of the modeled ionization ratios, though, are 
up to an order of magnitude lower than the Ulysses measure- 
ments, with slightly better agreement for the hotter Durham 
models. This discrepancy may be the result of our simple as- 
sumption of Maxwellian electron distributions and equal flow 
speeds for all ion charge states (see, e.g., Esser & Edgar 2000, 

2001) . A ubiquitous weak "halo" of suprathermal electrons in 
the low corona could be responsible for the larger overall de- 
gree of ionization that is observed. 

As one moves from the polar flux tube to the middle lat- 
itudes (i.e., from ps 750 km s -1 down to about 650 km 
s" 1 ) the modeled ionization ratio decreases, but then as one 
passes the bifurcation point between the low and high critical 
radii, the ionization ratio increases for the rest of the way from 
the mid-latitude region to the equator. The overall trend for 
slow wind to have a high ionization ratio is expected because 
of the higher temperatures in the narrow zone just above the 
transition region. In fact, the temperature at a constant height 
of z ~ O.O137?0 has the same quantitative trend as the mod- 
eled 7+ to 6+ ratios: there is an initial decrease from the 
polar model (0.37 MK) to just before the bifurcation latitude 
(0.28 MK), then an abrupt jump back to nearly the polar value 
(0.34 MK) followed by a steady increase toward the equato- 
rial flux tube (0.46 MK). The freezing in occurs at about this 
height, so it makes sense that the ionization state tracks these 
variations. Notice that this ionization ratio is sensitive only 
to the temperature trends at low heights (and that the 7+ to 
O fi+ ratio freezes in just above the transition region) and is not 
sensitive to the temperatures at larger heights in most of the 




FIG. 16. — FIP fractionation ratio of Fe to O abundance in units of their 
photospheric abundances. Symbols, line styles, and grayscale regions are the 
same as in Figure 15. 

extended corona. Some interpretations of the in situ charge 
states make the implicit assumption that high electron tem- 
peratures in the slow wind must persist all the way through 
the extended corona into interplanetary space, and our mod- 
els show that this need not be the case. 

2. A specific diagnostic of FIP fractionation has been 
computed for each of the standard and Durham model flux 
tubes. Figure 16 shows a commonly measured ratio of the 
abundances of a low-FIP element (Fe, 7.9 eV) and a moder- 
ately high-FIP element (O, 13.6 eV) normalized to their pho- 
tospheric abundance ratio. We solved equation d62i > as de- 
scribed in § 6, and by Laming (2004), for the various mod- 
els and compared the ratios to the Ulysses SWICS observa- 
tions of the Fe/O ratio. The SWICS ratio was normalized 
to an assumed photospheric iron-to-oxygen abundance ratio 
of 0.0468 (Grevesse & Sauval 1998) and the modeled quan- 
tity is naturally in units of the photospheric ratio. Because 
the FIP fractionation takes place in the uppermost part of the 
chromosphere (which is largely optically thin) and the very 
low transition region, we assumed "coronal" ionization equi- 
librium to hold for the modeled Fe and O ions. We used the 
code developed by Cranmer (2000) to compute the Fe and O 
ionization fractions that are used in equation d62l : the rates 
were taken largely from Mazzotta et al. (1998). 

Figure 16 shows that our standard pole-to-equator 
ZEPHYR model exhibits a roughly constant Fe/O ratio of 
about 2 for outflow speeds between 500 and 750 km s" 1 . This 
is in reasonable agreement with the Ulysses data, which have 
a roughly constant mean value of about 1 .5 over these speeds. 
The Durham model varies nonmonotonically between ratios 
of about 3 and 5.5 over its wider range of outflow speeds and 
is clearly inconsistent with most of the in situ data. The equa- 
torial Durham model, with ps 330 km s -1 , does seem to 
fall into agreement with the measured data, but the overall 
trend as a function of wind speed is not the same. The most 
promising models for matching the observed trend for a steep 
enhancement of the Fe/O ratio (as the wind speed decreases 
from 500 to 300 km s" 1 ) appear to be the active-region flux 
tubes discussed below in § 8.5. 

We should note that an alternate theoretical explanation for 
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FIP fractionation may also be consistent with the ZEPHYR 
models. Wang (1996) proposed that flux tubes with a stronger 
degree of downward heat conduction in the low corona could 
give rise to a larger amount of transient evaporative outflow 
of protons relative to neutral hydrogen atoms. These protons 
could collisionally dredge up singly ionized species (which 
are preferentially low-FIP) from the transition region and up- 
per chromosphere. In our models, the maximum value of 
iGcondl occurs within the transition region (at T « 10 5 K), and 
in the equatorial models this peak value is consistently larger 
than in the polar models. The ratio of equatorial to polar max- 
imum conduction rates is 1.7 for the standard model and 2.8 
for the Durham model. Whether this is enough to produce a 
FIP effect with Wang's (1996) mechanism is not known, but 
the variations are at least in the correct qualitative direction 
to produce more of an effect in the slow wind than in the fast 
wind. Note, though, that Schwadron et al. (1999) disagree 
with some of the basic assumptions of the Wang (1996) mech- 
anism, and a definitive explanation for the FIP effect remains 
elusive. 

There are additional observational diagnostics that appear 
to be correlated with fast and slow wind speeds. Kojima et 
al. (2004) and Suzuki (2006) found that the ratio of basal 
magnetic field strength to the superradial flux-tube divergence 
factor (i.e., Bbase//ss) may be a better predictor of solar wind 
speed than just using / ss itself (see, however, Wang & Shee- 
ley 2006). Mcintosh & Leamon (2005) found that the dif- 
ference in formation heights between the emission seen in 
the TRACE 1600 A and 1700 A filter bandpasses could be 
a useful diagnostic of the eventual wind speed in flux tubes 
that are associated with specific regions on the solar surface. 
Although there have been numerous observational inferences 
of blueshifts (i.e., outflow) in various spectral lines, it is un- 
clear to what degree these may be utilized as empirical "lower 
boundary conditions" for solar wind flux tubes (see reviews 
by Peter et al. 2004; Jones 2005). Finally, there are definite 
fast/slow wind correlations in the wealth of existing turbu- 
lence measurements — both in situ and remote sensing — that 
could be key tests of models that use waves and turbulence 
for energy and momentum deposition. Future development of 
the ZEPHYR models will encompass such additional compar- 
isons. 

8.5. Active Corona 

At times other than solar minimum, the coronal magnetic 
field is substantially more complex than the axisymmetric 
configuration modeled in the previous section. During the 
approach to solar maximum the large polar coronal holes dis- 
appear as magnetic flux of the opposite polarity is brought 
toward the poles. Active regions with strong fields emerge 
at all latitudes and give rise to an interconnected web of 
closed loops and intermittent open flux tubes. In the extended 
corona, active regions often have the appearance of cusp-like 
"active streamers" that seem like smaller, more intense ver- 
sions of the quiescent equatorial streamers that dominate at 
solar minimum (Newkirk 1967; Liewer et al. 2001; Ko et al. 
2002). There is increasing evidence that much of the slow 
solar wind is connected somehow to open flux tubes that are 
rooted in or near active regions (e.g., Neugebauer et al. 2002; 
Liewer et al. 2004). 

We investigated several different ways to model active- 
region flux tubes with the ZEPHYR code. First, we explored 
just changing the radial location of the cusp in the equato- 
rial Banaszkiewicz et al. (1998) flux tube. Observationally, 



streamers exhibit a wide range of cusp heights which gives 
rise to a large variation in their rates of superradial expansion. 
The cusp height was varied by first computing the ratio of the 
equatorial (6*o = 29.7°) to the polar (6a = 0) magnetic fields 
Bo(z), as shown in Figure 5, then shifting the height coordi- 
nate z either upward or downward by a constant multiplier. 
The "new" equatorial field was obtained by taking the prod- 
uct of this shifted ratio and the polar field strength. Changing 
the cusp height by up to a factor of two in either direction 
did not have much of an effect on the resulting solar wind 
properties. The mass loss rates varied by about 18% for these 
models (lower cusps having lower M), whereas the terminal 
speeds varied by less than 3%. This near constancy in the 
solar wind speed occurred despite the fact that / ss was made 
to increase and decrease by factors of ^2 about the equatorial 
value of 9.1. 

The above models cannot be considered true active- 
streamer models because they were not given the significantly 
stronger mean magnetic field strengths that active regions are 
observed to have in the lower atmosphere. On the smallest 
scales in the photosphere, plage and active regions exhibit a 
much higher filling factor of strong-field flux concentrations 
than quiet regions and coronal holes. In the most tightly- 
packed active regions these flux concentrations do not appear 
as isolated flux tubes but instead as fluted ribbons, buckled 
flux sheets, and flower-like patterns (e.g., Berger et al. 2004; 
Rouppe van der Voort et al. 2005). We are concerned here 
only with the subset of magnetic flux that extends out into in- 
terplanetary space, so we provisionally retain the "flux tube" 
picture. We thus modeled Ba(r) in an active-region flux tube 
by enhancing the field strength in the chromosphere and low 
corona. This is done in general agreement with the larger 
observed magnetic fields, but also to account for the smaller 
area A(r) subtended by open flux tubes that are "crowded" by 
neighboring closed loops in active regions. 

We added an active-field component to the standard equa- 
torial Banaszkiewicz et al. (1998) magnetic field by taking 



B (z) = max 



Be(z), 



(73) 



where Be(z) is the equatorial (8q = 29.7°) model and Ba is 
a constant that we set at 50 G. We produced a grid of mod- 
els that varied the scale height h between and 0.07 Rq. 
The model with /z — > is equivalent to the standard equa- 
torial model with no active-field enhancement. Equation 
d73l takes the greater of the two magnetic field quantities, 
which tends to retain the original Be in the photosphere, low 
chromosphere, and outer corona and wind, but produces the 
active-region enhancement in the upper chromosphere and 
low corona (see Figure 5). The smallest values of the scale 
height (h < 0.008 Rq) do not produce any enhancement in Bo 
because the quantity BA_e~ z ' h is rapidly decreasing and it never 
exceeds Be- For larger values of h, the magnetic field strength 
at the coronal base is enhanced and thus the Wang-Sheeley 
expansion factor / ss is enhanced as well. The maximum value 
of h = 0.07 R Q (i.e., 49 Mm) corresponds to / ss = 41. Even 
though we used the exact values of / ss computed for each 
model in the plots below, we found the following approximate 
fit to be a useful illustration of how / ss depends on h; 

9.1, Ii<0.02Rq 
9.1 + |3026(/j-0.02)| - 7 , h>0.02R & K ' 

where h is in units of Rq. The threshold value of 0.02 Rq 
corresponds to the point at which B^e" 2 -" 1 first produces a sig- 
nificant enhancement at Zbase = 0.04 Rq. 



/ss 



CORONAL HEATING AND SOLAR WIND ACCELERATION 



27 



1 

/ 

- 1 

i\ 

h 




>^ " ■ 


V 
!|l 




\. 


: j 


h - 


0.07 R Q ^ 




h = 


0.05 r q : 


ii I 


h - 


0.03 R Q ' 




.... h - 


0.02 R Q 


- 

l) 


h _ 


(std. equator) 


./ 

i i 


h - 


(std. polar) 
i i . .- 


0.01 


0.10 1.00 


10.00 100.00 




(r/R ) " 


- 1 



FIG . 17. — Radial dependence of temperature T for a selection of ZEPHYR 
models, highlighting the active-region extensions to the equatorial Bg(r) dis- 
cussed in § 8.5 (see captions). For comparison the standard-model polar (dot- 
ted line) and equatorial (dashed line) model temperatures are also shown. 

ZEPHYR models were computed for active region field en- 
hancements with h = 0.015, 0.02, 0.024, 0.03, 0.035, 0.04, 
0.045, 0.05, 0.06, and 0.07 R Q . The active region models ex- 
hibited lower terminal speeds and higher mass loss rates than 
the standard equatorial model. As h was increased from to 
0.07 Rq, Uoa decreased from 500 to about 390 km/s. The mass 
loss rate increased from 2.9 x 10" 14 to 3.5 x 10" 14 M Q yr _1 , 
thus increasing the proton number density at 1 AU from 6.4 
cm" 3 to about 10 cm" 3 . Figure 13 shows how decreases as 
/ ss increases for these models. There is reasonable agreement 
with the empirical Wang & Sheeley (1990) relationship, and 
it is important to reiterate that all models in Figure 13 were 
produced by fixing the basal parameters (v± & , L±q, F$q) at 
values appropriate for the polar coronal hole and varying only 
the magnetic field above the photosphere in order to produce 
the fast-to-slow wind transitions. 

Figures 15 and 16 illustrate the computed 7+ to O e+ 
freezing-in ratios and Fe/O (FIP fractionation) ratios for the 
active region models. These models mirror the observed 
trends for there to be larger ratios for both quantities at lower 
wind speeds. The strengthening of the FIP ratio, from ~ 1.5 to 
4 as the wind speed decreases from 500 to 400 km s" 1 , is espe- 
cially comparable to the observed increase toward the lowest 
speeds. 

Figure 17 compares the radial dependence of temperature 
for a selection of active models to that for the standard polar 
and equatorial models. As h increases from to 0.07 Rq, 
the peak coronal temperature increases from the equatorial 
model's value of 1.29 MK up to 1.67 MK. The heliocentric 
radius of this peak decreases from 1.26 to 1.11 Rq. These 
trends appear to be a continuation of the variations in coro- 
nal heating from the polar to the equatorial flux tubes. As 
one progresses from pole to equator (h = 0) to active region 
(h > 0), the models show a successively stronger decrease in 
the Alfven speed with height, more reflection (and damping) 
in the low corona, and thus less Alfven wave power available 
to heights at and above the critical point. 

The iterative convergence for the largest-/? active region 
ZEPHYR models was slower and less robust than for the mod- 



els discussed in §§ 8.1-8.4. The models having 0.035 < h < 
0.07 Rq were run for 400 outer iterations in order to reach ac- 
ceptable convergence parameter values: (SE) w 0.01 to 0.02. 
(Trial models with even larger values of h > 0.09 Rq did not 
converge below (SE) 0.04 after 400 iterations and were not 
used.) Unlike in Figure 4, where acceptable values of 0.01- 
0.02 were retained in nearly every iteration (after 100 or so 
outer iterations had passed), the active models exhibited these 
values only sporadically. Thus, the iteration with the low- 
est value of (SE) did not necessarily represent the "best" sin- 
gle solution. In order to evaluate the most likely time-steady 
plasma parameters, we considered a collection of the lowest- 
(SE) solutions sampled from the iteration stream for each 
model run. This was done by taking the distribution of values 
for a given parameter, such as u^, for all 400 outer iterations 
and weighting it by a factor proportional to (5E)~ 2 in order to 
deemphasize the early iterations that are farthest from conver- 
gence. Figure 13 shows both the mean and ±1 standard devi- 
ation values for this weighted distribution of well-converged 
terminal speed solutions. As h is increased, the standard devi- 
ation first decreases from about 20 to 6 km s" 1 , then increases 
again to 35 km s" 1 . 

It is unclear whether the nonconvergence of the large-/; ac- 
tive region models represents a failing of the numerical tech- 
nique or an intrinsic lack of a steady-state solution. The coro- 
nal sources of slow solar wind are traditionally viewed as 
more variable and filamentary than the source regions of fast 
wind. High-resolution coronagraphic observations of stream- 
ers reveal an almost continual release of low-contrast density 
inhomogeneities (Sheeley et al. 1997; Wang et al. 2000). Mul- 
tidimensional simulations of streamers have also exhibited 
time-variable solutions (e.g., Suess et al. 1996; Endeve et al. 
2003, 2004). For flux tubes in the vicinity of the boundary be- 
tween open and closed magnetic fields, the complex interplay 
between heat conduction and transverse pressure balance can 
lead to various instabilities that drive variability. Certain open 
magnetic field configurations (possibly like the one adopted 
in eq. 1731 ) may not be stable when considered in the context 
of its neighbors in a full-Sun model. 

Other techniques for constructing active-region solar wind 
models should be explored. The specific functional form 
of the Bo enhancement in the low corona can be improved 
in order to agree with, e.g., force-free magnetic field recon- 
structions of active regions (Schrijver et al. 2006). In ad- 
dition, the photospheric field strength may also be substan- 
tially higher in active regions than elsewhere. It is unclear at 
first glance how the photospheric wave fluxes would be mod- 
ified in these regions, though. For strong enough basal fields 
(i.e., in sunspots) convection is quenched and the wave fluxes 
should be lower. However, active regions exhibit a larger 
degree of magnetic flux emergence (e.g., Abramenko et al. 
2006) which could increase the frequency of impulsive hor- 
izontal "jumps" exhibited by flux-tube footpoints. As sum- 
marized above, flux tubes are more more highly concentrated 
in active-region intergranular lanes than elsewhere. The close 
presence of additional flux tubes may either increase the basal 
wave amplitude in a given tube (via "sympathetic" motions 
from the neighbors) or decrease it (because large horizontal 
displacements are blocked by the neighbors). Similarly, the 
correlation length L±q in these regions may be either lower 
(because of the smaller domain of motion when crowded) 
or higher (because the lane-filled "ribbons" may move to- 
gether). Comparisons between high-resolution observations 
of G-band bright points in active regions (Mostl et al. 2006) 
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and three-dimensional magnetoconvection simulations (e.g., 
Bushby & Houghton 2005) are needed in order to clarify these 
competing effects. 

9. SUMMARY AND DISCUSSION 

The primary aim of this paper has been to construct self- 
consistent models of chromospheric and coronal heating, 
waves and turbulence, and wind acceleration in an open mag- 
netic flux tube rooted in the solar photosphere. A key as- 
pect of the ZEPHYR models presented above is that the only 
true free parameters are: (1) the properties of the waves in- 
jected at the base, and (2) the background geometry and mag- 
netic field strength along the modeled flux tubes. Everything 
else (e.g., the radial dependence of the rates of chromospheric 
and coronal heating, the resulting temperature structure of the 
atmosphere, and the solar wind speed and mass flux) is an 
emergent property of the model and not an ad hoc input. In 
the above sense, the ZEPHYR models are similar to those of 
Suzuki & Inutsuka (2005, 2006); the differences between the 
two approaches involve mainly the specific physical processes 
that are assumed to dissipate the Alfven waves and heat the 
corona. 

As discussed in § 2, the ZEPHYR models presented here 
are limited by being one-dimensional, time-independent, and 
one-fluid solutions to the hydrodynamic conservation equa- 
tions, with heating that is derived only from a subset of all 
possible wave and turbulent dissipation processes. Even so, 
the results given in § 8 show that a realistic variation of 
asymptotic solar wind conditions can be produced by varying 
only the background magnetic field geometry, as predicted by 
Wang & Sheeley (1990, 1991, 2003, 2006). Specifically, our 
models show general agreement with some well known em- 
pirical correlations: i.e., a larger / ss (expansion factor) gives 
rise to a slower and denser wind, less intense Alfvenic fluctu- 
ations at 1 AU, and larger values of both the 7+ /0 6+ charge 
state ratio and the FlP-sensitive Fe/O abundance ratio. Satis- 
fying these kinds of scalings are necessary but not sufficient 
conditions for validating the idea that the wind is driven by a 
combination of MHD turbulence and non-WKB Alfven wave 
reflection. 

Future work must involve more physical realism for the 
models, expanded comparisons with existing observations, 
and predictions of as-yet unobserved quantities that may be 
key discriminators between competing theoretical models. An 
important component of all three efforts will be to model the 
divergent temperatures and flow speeds of protons, electrons, 
and various heavy ion species in the extended corona and he- 
liosphere. These mainly collisionless regions allow the ki- 
netic physics of wave dissipation to be probed to a level of 
detail not possible in a collisionally coupled (i.e., essentially 
one-fluid) plasma. Even in a perfectly collisional plasma, 
there can be macroscopic dynamical consequences depend- 
ing on how the energy is deposited into protons, electrons, 
and possibly heavy ions as well. For example, if all of the heat 
goes into electrons, there can be substantially more downward 
conduction than in a proton-heated model, which would affect 
the coronal temperature distribution (Hansteen & Leer 1995) 
and the stability of helmet streamers (Endeve et al. 2004). 
Two-fluid effects also may affect the phenomenological form 
of the MHD turbulent cascade assumed in the above models 
(e.g., Galtier 2006). 

The original motivation for this work was to understand 
the preferential ion heating and acceleration in the corona 
revealed by the UVCS (Ultraviolet Coronagraph Spectrom- 



eter) instrument on SOHO (Kohl et al. 1995, 1997, 2006). 
If MHD turbulence gives rise to cyclotron-resonant Alfven 
waves — which are likely to be responsible for the observed 
ion properties — then we must first have a large-scale de- 
scription of the energy flux injected into the turbulent cas- 
cade in order to further model the end-products of that cas- 
cade. We hope to link the macroscopic properties modeled 
by ZEPHYR with the microscopic kinetic processes studied 
by, e.g., Leamon et al. (1999, 2000), Cranmer & van Balle- 
gooijen (2003), Voitenko & Goossens (2003, 2004), Gary & 
Nishimura (2004), Dmitruk et al. (2004), Chandran (2005), 
Gary et al. (2006), and Markovskii et al. (2006). 

In addition to increasing the realism (and complexity) of 
the models, attempts should be made to extract the dominant 
physical processes of these models and create simpler "scal- 
ing laws" that can be used to estimate the wind conditions 
for varying magnetic field geometries and photospheric wave 
properties (see also Leer et al. 1982; Schwadron & McComas 
2003; Suzuki 2006). An obstacle to accomplishing this for the 
present ZEPHYR models is that the Alfven wave reflection is 
nontrivially coupled to the plasma state via the Alfven speed 
Va and solar wind speed u (and their gradients). The radial de- 
pendence of these quantities cannot be derived a priori from 
just the background magnetic field Bo(r). Approximations for 
Z_ and Z+ such as those given by Dmitruk et al. (2002) should 
be tested further to evaluate ways of estimating the amount 
of reflection without having to solve the full non-WKB equa- 
tions. 

We also intend to use the methodology developed here to 
model the winds of other kinds of stars on the cool side of the 
Hertzsprung-Russell (H-R) diagram. Similar ideas are begin- 
ning to be applied to younger and older stars with solar-type 
winds (e.g., Airapetian et al. 2000; Falceta-Goncalves et al. 
2006; Suzuki 2007; Holzwarth & Jardine 2007). A benefit of 
computing the Alfven wave evolution with a physically moti- 
vated damping rate (like eq. l47l ) is that the artificial "damp- 
ing lengths" used in the past are no longer needed. Many cool 
stars that exhibit observational evidence for winds have lower 
surface gravities than the Sun, stronger X-ray emission (i.e., 
more coronal heating), and much larger mass loss rates. These 
properties are likely to be causally linked, but no comprehen- 
sive and predictive models yet exist. 

In order to model cool-star winds with ZEPHYR we must 
have values for the input parameters that drive the atmo- 
spheric heating and wind acceleration. For example, we must 
be able to predict the properties of waves at the photospheric 
lower boundary for a given star. Amplitudes can be estimated 
from existing models of flux-tube wave generation from tur- 
bulent convection. Figure 18 shows an example of how trans- 
verse kink-mode wave amplitudes in stellar photospheres (for 
which we retain the symbol v^ ) can be estimated from the 
peak value of the subphotospheric convective velocity w Cjma x- 
Musielak et al. (2000) computed a series of stellar interior 
models with a range of effective temperatures (2000 to 10000 
K) and surface gravities (logg = 3, 4, and 5, with g in units 
of cm s~ 2 ) and presented the maximum values of the ratio 
Mc.max/ c s in their Figure 1 . Also, Musielak & Ulmschneider 
(2002) computed the photospheric transverse wave fluxes for 
the same grid of stellar models. We rescaled both quantities 
into velocity units using a sound speed consistent with T s s, a 
magnetic field strength consistent with the degree of pressure 
equipartition assumed by Musielak et al., and a photospheric 
density computed from the condition K^pH = 1, where H is 
the photospheric scale height (proportional to T e ff/g) and kr 
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FIG. 18. — Photospheric transverse velocity amplitude VX0 estimated 
from turbulent convection models of Musielak et al. (2000) and Musielak 
& Ulmschneider (2002), plotted against the maximum convective velocity 
Wc.max from below the photosphere. Models are for Population I stars with 
T c ff ranging from 2000 to 10000 K and logg = 5 (solid line), 4 (dashed line), 
and 3 (dotted line). The solar values used in this paper are denoted by the 
Sun symbol (0). 

is the Rosseland mean opacity interpolated from the same ta- 
bles used in other parts of the ZEPHYR code (see also Cran- 
mer 2005b). Note that although the two velocities span many 
orders of magnitude they "collapse" into something close to 
a one-to-one power-law relation (with vj_© cx k^), and that 
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the solar value used in § 8 coincides reasonably closely with 
the model curves. 

Another key input to cool-star wind models is an ade- 
quate description of the large-scale magnetic field. Observa- 
tional determinations of the field geometries of rapidly rotat- 
ing young stars have been made possible via Zeeman-Doppler 
imaging (e.g., Donati et al. 1990, 2003), time-resolved X-ray 
spectroscopy (Hussain et al. 2005), and combinations of these 
and related techniques. However, not enough is known about 
how the magnetic flux tubes break up — and penetrate the 
photosphere — on spatial scales relevant to convective gran- 
ulation and supergranulation. The relationship between the 
star-averaged mean field strength and the stronger fields in- 
ferred for thin flux tubes is not yet well understood (e.g., Saar 
2001). Ongoing improvements in the observations make pos- 
sible more comprehensive and predictive models. We antic- 
ipate that phenomenological approaches like those described 
in this paper can contribute to substantial progress in our un- 
derstanding of magnetic activity, coronal heating, and wind 
mass loss in a wide range of stars. 
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FIG. 6. — Contour plots of solar wind quantities resulting from varying the coronal heating parameters v ± q and L ± q\ (a) terminal wind speed Uoo in units of 
km s _1 , (b) mass loss rate M in units of Mq yr~' , (c) maximum coronal temperature in units of MK, (d) heliocentric critical radius in units of Rq . Also shown 
in each panel are the parameters chosen for the model of fast wind from a polar coronal hole discussed in § 8.3 (stars). 



